APIM - Control

Setup

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)

# Set options for analysis
shutdown = FALSE

options(dplyr.print_max = 100, 
        brms.backend = 'cmdstan',
        brms.file_refit = 'on_change' # "on_change" or "never"
        )

# functions for kable printing. 
print_df <- function(df, width = "auto") {
  kable(df) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "responsive"),
      full_width = F,
      position = "left",
      fixed_thead = T
    ) %>%
    column_spec(2:ncol(df), width = width) %>%
    row_spec(1:nrow(df), extra_css = "white-space: nowrap;") %>%
    scroll_box(width = '100%', height = '500px')
}

packing_ind <- function(df) {
  pack_rows(df,
    "Within-Person Effects", 2, 6) %>%
  pack_rows(
    "Between-Person Effects", 7, 10) %>%
  pack_rows(
    "Random Effects", 11, 13) %>%
  pack_rows("Additional Parameters", 14, 18)
}

packing_apim <- function(df) {
  pack_rows(df,
    "Within-Person Effects", 3, 12) %>%
  pack_rows(
    "Between-Person Effects", 13, 20) %>%
  pack_rows(
    "Random Effects", 21, 26) %>%
  pack_rows("Additional Parameters", 27, 31)
}


# Function to prepare for export
prepare_for_export <- function(model) {
  rep <- model
  rep$coefs <- rownames(rep)
  final <- cbind(rep$coefs, model)
  return(final)
}

Data Preparation

Scales and Scores

predictors_to_center <- c("persuasion",
                          'pressure',
                          'pushing', 
                          'weartime',
                          'barriers',
                          'support',
                          'comf',
                          'reas')


predictors_not_to_center <- c('day', 'userID', 'coupleID', # independent variables that don't need to be centered.
                              'aff', # outcomes (not centered)
                              'pa_sub',
                              'pa_obj', 
                              'anticipAff', 
                              'reactance', 
                              'plan',
                              'isWeekend',
                              'studyGroup',
                              'got_JITAI',
                              'skilled_support'
                              )



all_variables <- c(predictors_not_to_center, predictors_to_center)


df <- df %>% 
  mutate(ss_affect1 = ss_affect1 + 1, # re-code items from 0-5 to 1-6 scale.
         ss_affect2 = ss_affect2 + 1, 
         ss_affect3 = ss_affect3 + 1, 
         ss_affect4 = ss_affect4 + 1, 
         
         aff = (ss_affect1 + ss_affect3 + (7 - ss_affect2) + (7 - ss_affect4)) / 4, #Affect Scale

         
         pa_sub = ifelse(ss_pa == 0, 0, ss_pa_min_solo + ss_pa_min_collab),
         
         pa_obj = ifelse(wear_minutes > 600, minutes_mvpa_non_filtered, NA),
         
         day = day/54,
         
         barriers = (ss_barr_1 + ss_barr_2 + ss_barr_3 + ss_barr_4 + ss_barr_5 + ss_barr_6 + ss_barr_7) / 7,
      
         plan = ifelse(
           is.na(ss_pa_no), 
           ifelse(ss_pa_yes_0 == 1 | ss_pa_yes_1 == 1 | ss_pa_yes_2 == 1 | ss_pa_yes_3 == 1, 1,0),
           ifelse(ss_pa_no ==1, 1, 0)
         ), 
         
         got_JITAI = ifelse(((userID %% 2 != 0) & (trig_target_slot_0 %in% c("A", "C") | trig_target_slot_1 %in% c("A", "C") | trig_target_slot_2 %in% c("A", "C"))) |
                              ((userID %% 2 == 0) & (trig_target_slot_0 %in% c("B", "C") | trig_target_slot_1 %in% c("B", "C") | trig_target_slot_2 %in% c("B", "C"))), 1, 0),
         
         got_JITAI = ifelse(is.na(got_JITAI), 0, got_JITAI),
         
         skilled_support = ifelse(
           (studyGroup == 3 & ((day * 54 + 1) > 28)) | (studyGroup != 3 & ((day * 54 + 1) > 7)), 
           1, 0
           )
         
         )%>% 
  rename(persuasion = sp_psc_more,
         pressure = sp_nsc_more,
         pushing = sp_push_plan,
         anticipAff = ss_affect_con,
         reactance = ss_reactance,
         weartime = wear_minutes,
         # aff_reactance_persuasion = sp_persuasion_aff,
         #aff_reactance_pressure = sp_nsc_aff, # explorative
         #aff_reactance_pushing = sp_push_plan_aff, # explorative
         
         support = sp_emo_pleasure,
         comf = sp_emo_comf,
         reas = sp_emo_reass
         ) 

# pushing can only occur if there was a plan! Therefore, we set it to zero if there was no plan. 
df$pushing[df$plan == 0] <- 0

center Data

prepare_df <- function(df) {
  df_wide <- df %>%
    # center
    bmlm::isolate(by = "userID",
                  value = predictors_to_center,
                  which = "both") %>%
    # distinguish
    group_by(userID) %>%
    mutate(person_mean_pa_obj = mean(pa_obj, na.rm = TRUE)) %>%
    ungroup() %>%
    
    group_by(coupleID) %>%
    mutate(
      role = case_when(
        person_mean_pa_obj == max(person_mean_pa_obj, na.rm = TRUE) ~ 'A', 
        TRUE ~ 'B'
    )) %>%
    ungroup() %>%
  
    # Reshape the data to wide format
    pivot_wider(
      id_cols = c(coupleID, day),
      names_from = c(role),
      values_from = -c(1, 2, 4)
    )
    
    
  # Double Stacking Outcomes 
    
  df1 <- df_wide
  df2 <- df_wide
  
  df1$is_A <- 1
  df1$is_B <- 0
  df1$AorB <- 1
  df1$pa_sub <- df1$pa_sub_A
  df1$pa_obj <- df1$pa_obj_A
  df1$aff <- df1$aff_A
  df1$reactance <- df1$reactance_A
  df1$anticipAff <- df1$anticipAff_A
  #df1$aff_reactance_pressure <- df1$aff_reactance_pressure_A
  #df1$aff_reactance_pushing <- df1$aff_reactance_pushing_A
  
  df2$is_A <- 0
  df2$is_B <- 1
  df2$AorB <- 2
  df2$pa_sub <- df2$pa_sub_B
  df2$pa_obj <- df2$pa_obj_B
  df2$aff <- df2$aff_B
  df2$reactance <- df2$reactance_B
  df2$anticipAff <- df2$anticipAff_B
  #df2$aff_reactance_pressure <- df2$aff_reactance_pressure_B
  #df2$aff_reactance_pushing <- df2$aff_reactance_pushing_B
  
  
  df_double <- rbind(df1, df2)

  
  # Double Stacking predictors
  # persuasion
  df_double$persuasion_self <- df_double$persuasion_A * df_double$is_A + df_double$persuasion_B * df_double$is_B
  df_double$persuasion_partner <- df_double$persuasion_A * df_double$is_B + df_double$persuasion_B * df_double$is_A
  
  df_double$persuasion_self_cw <- df_double$persuasion_cw_A * df_double$is_A + df_double$persuasion_cw_B * df_double$is_B
  df_double$persuasion_partner_cw <- df_double$persuasion_cw_A * df_double$is_B + df_double$persuasion_cw_B * df_double$is_A
  
  df_double$persuasion_self_cb <- df_double$persuasion_cb_A * df_double$is_A + df_double$persuasion_cb_B * df_double$is_B
  df_double$persuasion_partner_cb <- df_double$persuasion_cb_A * df_double$is_B + df_double$persuasion_cb_B * df_double$is_A
  
  #pressure
  df_double$pressure_self <- df_double$pressure_A * df_double$is_A + df_double$pressure_B * df_double$is_B
  df_double$pressure_partner <- df_double$pressure_A * df_double$is_B + df_double$pressure_B * df_double$is_A
  
  df_double$pressure_self_cw <- df_double$pressure_cw_A * df_double$is_A + df_double$pressure_cw_B * df_double$is_B
  df_double$pressure_partner_cw <- df_double$pressure_cw_A * df_double$is_B + df_double$pressure_cw_B * df_double$is_A
  
  df_double$pressure_self_cb <- df_double$pressure_cb_A * df_double$is_A + df_double$pressure_cb_B * df_double$is_B
  df_double$pressure_partner_cb <- df_double$pressure_cb_A * df_double$is_B + df_double$pressure_cb_B * df_double$is_A
  
  # pushing
  df_double$pushing_self <- df_double$pushing_A * df_double$is_A + df_double$pushing_B * df_double$is_B
  df_double$pushing_partner <- df_double$pushing_A * df_double$is_B + df_double$pushing_B * df_double$is_A
  
  df_double$pushing_self_cw <- df_double$pushing_cw_A * df_double$is_A + df_double$pushing_cw_B * df_double$is_B
  df_double$pushing_partner_cw <- df_double$pushing_cw_A * df_double$is_B + df_double$pushing_cw_B * df_double$is_A
  
  df_double$pushing_self_cb <- df_double$pushing_cb_A * df_double$is_A + df_double$pushing_cb_B * df_double$is_B
  df_double$pushing_partner_cb <- df_double$pushing_cb_A * df_double$is_B + df_double$pushing_cb_B * df_double$is_A
  
  # Weartime
  
  df_double$weartime_self <- df_double$weartime_A * df_double$is_A + df_double$weartime_B * df_double$is_B
  df_double$weartime_partner <- df_double$weartime_A * df_double$is_B + df_double$weartime_B * df_double$is_A
  
  df_double$weartime_self_cw <- df_double$weartime_cw_A * df_double$is_A + df_double$weartime_cw_B * df_double$is_B
  df_double$weartime_partner_cw <- df_double$weartime_cw_A * df_double$is_B + df_double$weartime_cw_B * df_double$is_A
  
  df_double$weartime_self_cb <- df_double$weartime_cb_A * df_double$is_A + df_double$weartime_cb_B * df_double$is_B
  df_double$weartime_partner_cb <- df_double$weartime_cb_A * df_double$is_B + df_double$weartime_cb_B * df_double$is_A
  
  # barriersiers
  df_double$barriers_self <- df_double$barriers_A * df_double$is_A + df_double$barriers_B * df_double$is_B
  df_double$barriers_partner <- df_double$barriers_A * df_double$is_B + df_double$barriers_B * df_double$is_A
  
  df_double$barriers_self_cw <- df_double$barriers_cw_A * df_double$is_A + df_double$barriers_cw_B * df_double$is_B
  df_double$barriers_partner_cw <- df_double$barriers_cw_A * df_double$is_B + df_double$barriers_cw_B * df_double$is_A
  
  df_double$barriers_self_cb <- df_double$barriers_cb_A * df_double$is_A + df_double$barriers_cb_B * df_double$is_B
  df_double$barriers_partner_cb <- df_double$barriers_cb_A * df_double$is_B + df_double$barriers_cb_B * df_double$is_A
  
  
  
  ## For Sensitivity Analyses necessary
  
  # Plans
  df_double$plan_self <- as.numeric(df_double$plan_A) * df_double$is_A + as.numeric(df_double$plan_B) * df_double$is_B
  df_double$plan_partner <- as.numeric(df_double$plan_A) * df_double$is_B + as.numeric(df_double$plan_B) * df_double$is_A
  
  # Support
  df_double$support_self <- df_double$support_A * df_double$is_A + df_double$support_B * df_double$is_B
  df_double$support_partner <- df_double$support_A * df_double$is_B + df_double$support_B * df_double$is_A
  
  df_double$support_self_cw <- df_double$support_cw_A * df_double$is_A + df_double$support_cw_B * df_double$is_B
  df_double$support_partner_cw <- df_double$support_cw_A * df_double$is_B + df_double$support_cw_B * df_double$is_A
  
  df_double$support_self_cb <- df_double$support_cb_A * df_double$is_A + df_double$support_cb_B * df_double$is_B
  df_double$support_partner_cb <- df_double$support_cb_A * df_double$is_B + df_double$support_cb_B * df_double$is_A
  
  # JITAIs
  df_double$got_JITAI_self <- df_double$got_JITAI_A * df_double$is_A + df_double$got_JITAI_B * df_double$is_B
  df_double$got_JITAI_partner <- df_double$got_JITAI_A * df_double$is_B + df_double$got_JITAI_B * df_double$is_A  
  
  
  df_double$pa_obj_log <- log(df_double$pa_obj)
      
  return(list(df_wide, df_double)) 
}


# And prepare original dataset
df_list <- prepare_df(df)

df_wide <- df_list[[1]]
df_double <- df_list[[2]]

Creating function to report measures

report_measures <- function(data, measures, ICC = TRUE) {
  
  if (ICC & any(sapply(data[,measures], is.factor))) {
    warning("Cannot add ICCs when Factors are in Dataframe.")
    ICC <- FALSE
  }
  
  measures_table <- as.data.frame(report::report_table(data[, measures]))

  measures_table <- measures_table %>%
    mutate(across(.cols = c('Mean', 'SD'),
                  .fns = ~format(round(.x, 2), nsmall = 2)))
  
  
  measures_table$Missing <- ifelse(
    is.na(measures_table$percentage_Missing),
    NA,
    paste0(round(measures_table$percentage_Missing), '%')
  )
  
  measures_table$Range <- ifelse(
    is.na(measures_table$Min) | is.na(measures_table$Max),
    NA,
    paste0(format(round(measures_table$Min, 2), nsmall = 2), '-', format(round(measures_table$Max, 2), nsmall = 2))
  )
    
  
  
  finished_df <- measures_table
  
  if (ICC) {
    cors <- wbCorr(data[,measures], data$coupleID)
    
    finished_df <- cbind(
      measures_table[, c('Variable', 'n_Obs', 'Missing', 'Mean', 'SD', 'Range')],
      ICC = format(round(get_icc(cors)$ICC, 2), nsmall = 2)
    )
  } else {
    measures_table$percentage_Obs <- ifelse(
      is.na(measures_table$percentage_Obs),
      NA,
      paste0(round(measures_table$percentage_Obs), '%')
    )
    
    finished_df <- measures_table[, c('Variable', 'Level', 'n_Obs', 'percentage_Obs', 'Missing', 'Mean', 'SD', 'Range')]
  }
  
  # Make them all Text for easier copying to excel. 
  finished_df[] <- lapply(finished_df, as.character)
  
  return(finished_df)
}

Sample Description

# mean of both relationship durations reported
df_wide$rel_duration_A <- df_wide$pre_rel_duration_y_A + df_wide$pre_rel_duration_m_A / 12 # in months
df_wide$rel_duration_B <- df_wide$pre_rel_duration_y_B + df_wide$pre_rel_duration_m_B / 12
df_wide$Relationship_Duration <- (df_wide$rel_duration_A + df_wide$rel_duration_B) / 2

# mean of both cohabiting duration reported
df_wide$hab_duration_A <- df_wide$pre_hab_duration_y_A + df_wide$pre_hab_duration_m_A / 12
df_wide$hab_duration_B <- df_wide$pre_hab_duration_y_B + df_wide$pre_hab_duration_m_B / 12 
df_wide$Cohabiting_Duration <- (df_wide$hab_duration_A + df_wide$hab_duration_B) / 2

#Gender of Each Partner
df_wide$Gender_A <- factor(df_wide$gender_A, levels = c(1,2,3), labels = c('Male','Female', 'Other'))
df_wide$Gender_B <- factor(df_wide$gender_B, levels = c(1,2,3), labels = c('Male','Female', 'Other'))
#Same sex couples or not
df_wide$Same_Sex <- factor(df_wide$Gender_A == df_wide$Gender_B, levels = c(FALSE, TRUE), labels = c('Same-Sex Couple', 'Mixed-Sex Couple'))

#Handedness
df_wide$Handedness_A <- factor(df_wide$pre_handedness_A, levels = c(0, 1, 2), c('Right','Left', 'Ambidextrous'))
df_wide$Handedness_B <- factor(df_wide$pre_handedness_B, levels = c(0, 1, 2), c('Right','Left', 'Ambidextrous'))



# Income (mean of both partner's report)

merge_income <- function(income1, income2) {
  merged_income <- numeric(length(income1))
  
  # Loop through each pair of incomes
  for (i in seq_along(income1)) {
    # Handle NA
    if (is.na(income1[i])) {
      merged_income[i] <- income2[i]
    } 
    
    else if (is.na(income2[i])) {
      merged_income[i] <- income1[i]
    }
    
    # if both are informative, take mean and round
    else if (income1[i] %in% 1:6 && income2[i] %in% 1:6) {
      merged_income[i] <- round((income1[i] + income2[i]) / 2)
    }
    
    # if one is informative and the other not, use the informative one. 
    else if (income1[i] %in% 1:6) {
      merged_income[i] <- income1[i]
    }
    else if (income2[i] %in% 1:6) {
      merged_income[i] <- income2[i]
    }
    
    # Now we only have cases, where both are either 70 or 99. We simply report "undisclosed". 
    else {
      merged_income[i] <- 99
    }
  }
  
  return(merged_income)
}


# Apply the function
numeric_income <- merge_income(df_wide$pre_income_1_A, df_wide$pre_income_1_B)

# Convert to factor
df_wide$Household_Income <- factor(numeric_income, levels = c(1,2,3,4,5,6,70,99), labels = c(
    "up to CHF 2'000.-", 
    "CHF 2'001.- to CHF 4'000.-",
    "CHF 4'001.- to CHF 6'000.-",
    "CHF 6'001.- to CHF 8'000.-",
    "CHF 8'001.- to CHF 10'000.-", 
    "above CHF 10'000.-", 
    "I don't know",
    "Undisclosed"
))





# Education for both separate
df_wide$Highest_Education_A <- factor(df_wide$pre_education_A, levels = c(1,2,3,4,5,6,7), labels = c(
  "(still) no school diploma",
  "compulsory education (9 years)",
  "vocational training (apprenticeship)",
  "Matura (university entrance qualification)",
  "Bachelor's degree", 
  "Master's degree",
  "Doctorate degree"
  )
)

df_wide$Highest_Education_B <- factor(df_wide$pre_education_B, levels = c(1,2,3,4,5,6,7), labels = c(
  "(still) no school diploma",
  "compulsory education (9 years)",
  "vocational training (apprenticeship)",
  "Matura (university entrance qualification)",
  "Bachelor's degree", 
  "Master's degree",
  "Doctorate degree"
  )
)



# Marital Status (married)
df_wide$Marital_Status_A <- df_wide$pre_mat_stat_A == 1
df_wide$Marital_status_B <- df_wide$pre_mat_stat_B == 1

# if at least one of them said they are married, we go for married
df_wide$Marital_Status <- factor((df_wide$Marital_Status_A + df_wide$Marital_status_B) > 0, levels = c(FALSE, TRUE), labels = c('Not Married', 'Married'))



# Children
df_wide$Children <- factor((df_wide$pre_child_option_A + df_wide$pre_child_option_B) > 0, levels = c(FALSE, TRUE), labels = c('Have Children', 'No Children'))
# number of chilren (ONLY FOR COUPLES THAT HAVE CHILDREN)
df_wide$nChildren_couples_with_children <- (df_wide$pre_child_nr_childs_A + df_wide$pre_child_nr_childs_B) / 2 
# number of children if every couple is considered. 
df_wide$nChildren_all_couples <- df_wide$nChildren_couples_with_children
df_wide$nChildren_all_couples[is.na(df_wide$nChildren_all_couples)] <- 0


# BMI (kg/m^2)
df_wide$pre_height_A <- df_wide$pre_height_A / 100 # to meters
df_wide$pre_height_B <- df_wide$pre_height_B / 100 # to meters

df_wide$BMI_A <- df_wide$pre_weight_A / (df_wide$pre_height_A^2)
df_wide$BMI_B <- df_wide$pre_weight_B / (df_wide$pre_height_B^2)



sample_measures <- c(
  "Gender_A", "Gender_B", "Same_Sex",
  "pre_age_A", "pre_age_B", "BMI_A", "BMI_B", "Handedness_A", "Handedness_B", # physical stats
  "Household_Income", "Highest_Education_A", "Highest_Education_B", # SES
  "Marital_Status", "Relationship_Duration", "Cohabiting_Duration",
  "Children", "nChildren_all_couples", "nChildren_couples_with_children" 
  )



sample_table <- report_measures(df_wide, sample_measures, ICC=F)


rownames(sample_table) <- NULL

openxlsx::write.xlsx(sample_table, 'Output/AsPreregistred/Sample.xlsx')

print_df(sample_table)

Variable	Level	n_Obs	percentage_Obs	Missing	Mean	SD	Range
Gender_A	Other	0	0%	NA	NA	NA	NA
Gender_A	Female	825	39%	NA	NA	NA	NA
Gender_A	Male	1265	61%	NA	NA	NA	NA
Gender_B	Other	0	0%	NA	NA	NA	NA
Gender_B	Male	825	39%	NA	NA	NA	NA
Gender_B	Female	1265	61%	NA	NA	NA	NA
Same_Sex	Mixed-Sex Couple	110	5%	NA	NA	NA	NA
Same_Sex	Same-Sex Couple	1980	95%	NA	NA	NA	NA
pre_age_A	NA	2090	NA	0%	34.05	11.19	19.00-58.00
pre_age_B	NA	2090	NA	0%	33.97	10.72	19.00-60.00
BMI_A	NA	2090	NA	0%	24.58	3.77	16.37-33.95
BMI_B	NA	2090	NA	0%	25.29	4.33	17.32-33.95
Handedness_A	Ambidextrous	55	3%	NA	NA	NA	NA
Handedness_A	Left	330	16%	NA	NA	NA	NA
Handedness_A	Right	1705	82%	NA	NA	NA	NA
Handedness_B	Ambidextrous	0	0%	NA	NA	NA	NA
Handedness_B	Left	275	13%	NA	NA	NA	NA
Handedness_B	Right	1815	87%	NA	NA	NA	NA
Household_Income	I don’t know	0	0%	NA	NA	NA	NA
Household_Income	up to CHF 2’000.-	110	5%	NA	NA	NA	NA
Household_Income	CHF 2’001.- to CHF 4’000.-	165	8%	NA	NA	NA	NA
Household_Income	CHF 6’001.- to CHF 8’000.-	165	8%	NA	NA	NA	NA
Household_Income	CHF 4’001.- to CHF 6’000.-	275	13%	NA	NA	NA	NA
Household_Income	Undisclosed	275	13%	NA	NA	NA	NA
Household_Income	CHF 8’001.- to CHF 10’000.-	440	21%	NA	NA	NA	NA
Household_Income	above CHF 10’000.-	660	32%	NA	NA	NA	NA
Highest_Education_A	(still) no school diploma	0	0%	NA	NA	NA	NA
Highest_Education_A	compulsory education (9 years)	0	0%	NA	NA	NA	NA
Highest_Education_A	Doctorate degree	0	0%	NA	NA	NA	NA
Highest_Education_A	vocational training (apprenticeship)	330	16%	NA	NA	NA	NA
Highest_Education_A	Master’s degree	495	24%	NA	NA	NA	NA
Highest_Education_A	Matura (university entrance qualification)	605	29%	NA	NA	NA	NA
Highest_Education_A	Bachelor’s degree	660	32%	NA	NA	NA	NA
Highest_Education_B	(still) no school diploma	0	0%	NA	NA	NA	NA
Highest_Education_B	compulsory education (9 years)	0	0%	NA	NA	NA	NA
Highest_Education_B	Doctorate degree	0	0%	NA	NA	NA	NA
Highest_Education_B	vocational training (apprenticeship)	110	5%	NA	NA	NA	NA
Highest_Education_B	Bachelor’s degree	605	29%	NA	NA	NA	NA
Highest_Education_B	Master’s degree	660	32%	NA	NA	NA	NA
Highest_Education_B	Matura (university entrance qualification)	715	34%	NA	NA	NA	NA
Marital_Status	Married	715	34%	NA	NA	NA	NA
Marital_Status	Not Married	1375	66%	NA	NA	NA	NA
Relationship_Duration	NA	2090	NA	0%	9.23	9.03	0.58-36.00
Cohabiting_Duration	NA	2090	NA	0%	7.53	9.14	0.25-33.00
Children	No Children	605	29%	NA	NA	NA	NA
Children	Have Children	1485	71%	NA	NA	NA	NA
nChildren_all_couples	NA	2090	NA	0%	0.47	0.85	0.00- 3.00
nChildren_couples_with_children	NA	2090	NA	74%	1.80	0.60	1.00- 3.00

Descriptives of Measures

Main Measures

# Main Measures
main_constructs <- c("persuasion_A", "persuasion_B", 
                     "pressure_A", "pressure_B", 
                     "pushing_A", "pushing_B",
                     
                     "aff_A", "aff_B", 
                     "pa_sub_A", "pa_sub_B","pa_obj_A", "pa_obj_B", 
                     "anticipAff_A", "anticipAff_B", 
                     "reactance_A", "reactance_B"
                     )

main_descriptives <- report_measures(df_wide, main_constructs)

openxlsx::write.xlsx(main_descriptives, 'Output/AsPreregistred/DescriptivesMain.xlsx')

print_df(main_descriptives)

	Variable	n_Obs	Missing	Mean	SD	Range	ICC
6	persuasion_A	2090	7%	0.43	1.09	0.00- 5.00	0.24
5	persuasion_B	2090	5%	0.41	1.08	0.00- 5.00	0.23
2	pressure_A	2090	7%	0.13	0.62	0.00- 5.00	0.63
1	pressure_B	2090	5%	0.12	0.62	0.00- 5.00	0.53
4	pushing_A	2090	7%	0.18	0.68	0.00- 5.00	0.09
3	pushing_B	2090	5%	0.15	0.64	0.00- 5.00	0.12
12	aff_A	2090	6%	4.89	1.13	1.00- 6.00	0.47
11	aff_B	2090	5%	4.77	1.14	1.00- 6.00	0.36
14	pa_sub_A	2090	6%	31.66	56.73	0.00-720.00	0.17
13	pa_sub_B	2090	5%	29.16	52.77	0.00-480.00	0.15
16	pa_obj_A	2090	12%	164.63	133.55	7.75-971.25	0.14
15	pa_obj_B	2090	11%	124.32	95.63	5.75-925.25	0.20
10	anticipAff_A	2090	6%	1.50	2.23	-5.00- 5.00	0.39
9	anticipAff_B	2090	5%	1.25	2.39	-5.00- 5.00	0.37
7	reactance_A	2090	80%	0.70	1.24	0.00- 4.00	0.52
8	reactance_B	2090	82%	0.90	1.39	0.00- 5.00	0.42

Including all Covariates

# With Covariates


# JITAIs as factor
df_wide$got_JITAI_A_factor <- factor(df_wide$got_JITAI_A, levels = c(0,1), labels = c("No JITAI", "Got a JITAI"))
df_wide$got_JITAI_B_factor <- factor(df_wide$got_JITAI_B, levels = c(0,1), labels = c("No JITAI", "Got a JITAI"))


# weekend as factor
df_wide$isWeekend_factor <- factor(df_wide$isWeekend_A, levels = c(0,1), labels = c("weekday", "weekend day"))


# Plan as a factor
df_wide$plan_A_factor <- factor(df_wide$plan_A, levels = c(0,1), labels = c("No Plan", "Yes Plan"))
df_wide$plan_B_factor <- factor(df_wide$plan_B, levels = c(0,1), labels = c("No Plan", "Yes Plan"))


# Study Group as a Factor
df_wide$studyGroup_factor <- factor(df_wide$studyGroup_A - 1, levels = c(0,1,2), labels = c("Allways Interventions", "First 3 weeks", "last 3 weeks"))

# Skilled support as factor
df_wide$skilled_support_factor <- factor(df_wide$skilled_support_A, levels = c(0,1), labels = c("Did not yet have SS", "Did already have SS"))



all_constructs <- c(
  main_constructs,
  "day",
  "weartime_A", "weartime_B",
  "barriers_A" , "barriers_B",
  "isWeekend_factor",
  "plan_A_factor", "plan_B_factor",
  "studyGroup_factor",
  "support_A", "support_B",
  "got_JITAI_A_factor", "got_JITAI_B_factor",
  "skilled_support_factor"
)


all_descriptives <- report_measures(df_wide, all_constructs, ICC = F)

openxlsx::write.xlsx(all_descriptives, 'Output/AsPreregistred/DescriptivesAll.xlsx')

print_df(all_descriptives)

	Variable	Level	n_Obs	percentage_Obs	Missing	Mean	SD	Range
29	persuasion_A	NA	2090	NA	7%	0.43	1.09	0.00- 5.00
30	persuasion_B	NA	2090	NA	5%	0.41	1.08	0.00- 5.00
31	pressure_A	NA	2090	NA	7%	0.13	0.62	0.00- 5.00
32	pressure_B	NA	2090	NA	5%	0.12	0.62	0.00- 5.00
33	pushing_A	NA	2090	NA	7%	0.18	0.68	0.00- 5.00
34	pushing_B	NA	2090	NA	5%	0.15	0.64	0.00- 5.00
18	aff_A	NA	2090	NA	6%	4.89	1.13	1.00- 6.00
19	aff_B	NA	2090	NA	5%	4.77	1.14	1.00- 6.00
27	pa_sub_A	NA	2090	NA	6%	31.66	56.73	0.00- 720.00
28	pa_sub_B	NA	2090	NA	5%	29.16	52.77	0.00- 480.00
25	pa_obj_A	NA	2090	NA	12%	164.63	133.55	7.75- 971.25
26	pa_obj_B	NA	2090	NA	11%	124.32	95.63	5.75- 925.25
20	anticipAff_A	NA	2090	NA	6%	1.50	2.23	-5.00- 5.00
21	anticipAff_B	NA	2090	NA	5%	1.25	2.39	-5.00- 5.00
35	reactance_A	NA	2090	NA	80%	0.70	1.24	0.00- 4.00
36	reactance_B	NA	2090	NA	82%	0.90	1.39	0.00- 5.00
24	day	NA	2090	NA	0%	0.50	0.29	0.00- 1.00
39	weartime_A	NA	2090	NA	2%	1036.85	374.95	0.00-1440.00
40	weartime_B	NA	2090	NA	0%	1077.63	392.29	0.00-1440.00
22	barriers_A	NA	2090	NA	54%	0.04	1.58	-3.29- 5.00
23	barriers_B	NA	2090	NA	57%	0.10	1.70	-3.57- 5.00
6	isWeekend_factor	weekend day	608	29%	NA	NA	NA	NA
14	isWeekend_factor	weekday	1482	71%	NA	NA	NA	NA
2	plan_A_factor	missing	136	7%	NA	NA	NA	NA
11	plan_A_factor	Yes Plan	959	46%	NA	NA	NA	NA
12	plan_A_factor	No Plan	995	48%	NA	NA	NA	NA
1	plan_B_factor	missing	102	5%	NA	NA	NA	NA
10	plan_B_factor	Yes Plan	901	43%	NA	NA	NA	NA
13	plan_B_factor	No Plan	1087	52%	NA	NA	NA	NA
7	studyGroup_factor	last 3 weeks	660	32%	NA	NA	NA	NA
8	studyGroup_factor	Allways Interventions	715	34%	NA	NA	NA	NA
9	studyGroup_factor	First 3 weeks	715	34%	NA	NA	NA	NA
37	support_A	NA	2090	NA	7%	0.87	1.45	0.00- 5.00
38	support_B	NA	2090	NA	5%	0.94	1.58	0.00- 5.00
3	got_JITAI_A_factor	Got a JITAI	287	14%	NA	NA	NA	NA
17	got_JITAI_A_factor	No JITAI	1803	86%	NA	NA	NA	NA
4	got_JITAI_B_factor	Got a JITAI	298	14%	NA	NA	NA	NA
16	got_JITAI_B_factor	No JITAI	1792	86%	NA	NA	NA	NA
5	skilled_support_factor	Did not yet have SS	518	25%	NA	NA	NA	NA
15	skilled_support_factor	Did already have SS	1572	75%	NA	NA	NA	NA

Correlations Main Constructs

cors <- wbCorr(df_wide[,c(main_constructs)], df_wide$coupleID, method = 'spearman')

main_cors <- summary(cors, 'wb')$merged_wb


openxlsx::write.xlsx(main_cors, 'Output/AsPreregistred/Correlations.xlsx')

print_df(main_cors, width = '7em')

	persuasion_A	persuasion_B	pressure_A	pressure_B	pushing_A	pushing_B	aff_A	aff_B	pa_sub_A	pa_sub_B	pa_obj_A	pa_obj_B	anticipAff_A	anticipAff_B	reactance_A	reactance_B
persuasion_A	[0.24]	0.23***	0.29***	0.19***	0.38***	0.07**	-0.02	0.03	0.14***	0.11***	0.01	0.00	0.02	0.03	-0.08	-0.01
persuasion_B	0.32	[0.23]	0.03	0.27***	0.08***	0.37***	-0.01	-0.01	0.20***	0.20***	0.07**	0.13***	0.05*	0.04	0.01	0.02
pressure_A	0.46**	0.02	[0.63]	0.30***	0.33***	-0.00	0.01	0.01	-0.01	-0.01	-0.03	-0.01	0.07**	0.06**	0.24***	0.03
pressure_B	0.47**	0.34*	0.45**	[0.53]	0.05*	0.23***	-0.02	-0.02	0.02	-0.04*	-0.08**	0.06*	0.04	0.02	0.05	0.28***
pushing_A	0.55***	0.23	0.48**	0.24	[0.09]	0.17***	-0.01	0.02	0.10***	0.10***	0.02	-0.02	0.02	0.02	0.09	-0.11
pushing_B	0.24	0.67***	0.11	0.38*	0.35*	[0.12]	0.03	0.02	0.17***	0.13***	0.07**	0.10***	0.04	0.04	0.01	0.00
aff_A	0.35*	0.19	-0.03	-0.06	0.24	-0.03	[0.47]	0.19***	0.19***	0.11***	0.13***	0.12***	0.30***	0.10***	-0.02	0.00
aff_B	0.23	-0.04	0.01	-0.12	0.35*	-0.01	0.41*	[0.36]	0.10***	0.20***	0.10***	0.18***	0.08***	0.36***	-0.02	-0.04
pa_sub_A	0.12	0.23	-0.27	-0.06	0.13	0.20	0.54***	0.47**	[0.17]	0.56***	0.28***	0.24***	0.19***	0.11***	-0.09	0.02
pa_sub_B	0.08	0.26	-0.16	-0.19	0.13	0.11	0.19	0.46**	0.70***	[0.15]	0.22***	0.33***	0.11***	0.20***	0.04	-0.05
pa_obj_A	0.03	0.22	0.02	-0.22	0.26	0.20	0.12	0.25	0.30	0.38*	[0.14]	0.23***	0.05*	0.03	0.01	0.03
pa_obj_B	-0.06	0.18	0.04	-0.20	0.19	0.10	-0.00	0.40*	0.37*	0.54***	0.76***	[0.20]	0.11***	0.11***	0.11*	0.13*
anticipAff_A	0.24	0.40*	-0.35*	-0.06	0.14	0.20	0.47**	0.38*	0.66***	0.49**	0.40*	0.39*	[0.39]	0.15***	-0.05	0.04
anticipAff_B	0.44**	0.12	0.07	0.25	0.39*	0.17	0.10	0.62***	0.29	0.44**	0.11	0.28	0.41**	[0.37]	0.00	-0.01
reactance_A	0.07	0.03	0.32*	0.07	0.14	0.03	-0.16	0.03	-0.13	-0.09	0.40*	0.37*	0.03	0.04	[0.52]	0.26***
reactance_B	0.17	0.30	-0.01	0.22	-0.16	0.29	-0.04	-0.15	-0.01	-0.03	0.07	0.06	0.11	0.02	0.01	[0.42]

plot(cors,'w')

## This may take a while...

#Modelling

Writing custion summary functions to facilitate reporting

summarize_brms <- function(model, short_version = FALSE, exponentiate = FALSE, model_rows_fixed = NULL, model_rows_random = NULL) {

  summ_og <- summary(model)
  summ <- summ_og$fixed
  rands <- summ_og$random[[1]][grep('sd\\(', rownames(summ_og$random[[1]])), ]
  
  ar1sterr <- summ_og$cor_pars
  sdresid <- summ_og$spec_pars
  
  rands <- rbind(rands, ar1sterr, sdresid)
  
  # Get rows that we need
  if (is.null(model_rows_fixed)) {
    model_rows_fixed <- rownames(summ)
  }
  if (is.null(model_rows_random)) {
    model_rows_random <- rownames(rands)
  }
  
  summ <- summ[model_rows_fixed, ]
  rands <- format(round(rands[model_rows_random, ], 2), nsmall = 2)
  
  Rhat <- format(round(summ$Rhat, 3), nsmall = 3)

  summ$Est.Error <- format(round(summ$Est.Error, 2), nsmall=2)
  summ$Bulk_ESS <- format(round(summ$Bulk_ESS, 2), nsmall=2)
  summ$Tail_ESS <- format(round(summ$Tail_ESS, 2), nsmall=2)
  
  #exponentiate if needed
  if (exponentiate) {
    summ$Estimate <- exp(summ$Estimate)
  }
    
  # Add stars based on Credible Interval
  significance <- ifelse(
    (summ$`l-95% CI` > 0 & summ$`u-95% CI` > 0) | 
    (summ$`l-95% CI` < 0 & summ$`u-95% CI` < 0), 
    '*', 
    '')
  
  # exponentiate CI
  if (exponentiate) {
    summ$`u-95% CI`<- exp(summ$`u-95% CI`)
    summ$`l-95% CI` <- exp(summ$`l-95% CI`)
  }
      

  
  # Round CI
  summ$`l-95% CI` <- format(round(summ$`l-95% CI`, 2), nsmall = 2)
  summ$`u-95% CI` <- format(round(summ$`u-95% CI`, 2), nsmall = 2)
  
  
  summ$Estimate <- ifelse( 
    is.na(summ$Estimate), 
    NA, 
    paste0(
      format(round(summ$Estimate,2), nsmall=2), 
      significance
    )
  )
  
  
  # Rename "Estimate" column correctly
  if (exponentiate) {
    if ('bernoulli' %in% model$family) {
      correct_name <- 'OR'
    } else if ('negbinomial' %in% model$family) {
      correct_name <- 'IRR'
    } else {
      correct_name <- 'exp(Est.)'
      warning(
        "Coefficients were exponentiated. Double check if this was intended."
        )
    }
    
  } else {
    correct_name <- 'b'
  }
  
  colnames(summ)[1] <- correct_name
  colnames(rands)[1] <- correct_name
  
  
  # Add rhat rounded to 3 digits back
  summ$Rhat <- Rhat
  
  # Add back all rownames
  rownames(summ) <- model_rows_fixed
  rownames(rands) <- model_rows_random
  
  # Add random effects to fixed effects df. 
  fullrep <- rbind(summ, rands)
  

  
  # We are finished, if we want the full table
  if (!short_version) {
      return(fullrep[ , c(1, 3, 4, 5, 6, 7)])
  }
  
  
  fullrep$`95% CI` <-  ifelse(
    is.na(fullrep[correct_name]) | fullrep$`u-95% CI` == "NA",
    NA,
    paste0("[", fullrep$`l-95% CI`, ", ", fullrep$`u-95% CI`, "]")
  )

  return(fullrep[ , c(1, 8)])
}


report_side_by_side <- function(..., model_rows_random, model_rows_fixed) {
  models <- list(...)
  model_names <- sapply(substitute(list(...))[-1], deparse)
  names(models) <- model_names

  side_by_side <- NULL
  for (i in seq_along(models)) {
    model <- models[[i]]
    model_name <- model_names[i]
    print(model_name)
    if ('bernoulli' %in% model$family | 'negbinomial' %in% model$family | grepl('log', model_name)) {
      exponentiate <- TRUE
    } else {
      exponentiate <- FALSE
    }
    model_summary <- summarize_brms(
      model, short_version = TRUE, 
      exponentiate = exponentiate, 
      model_rows_random = model_rows_random,
      model_rows_fixed = model_rows_fixed
      )
    
    colnames(model_summary) <- paste(colnames(model_summary), model_name)
    
    if (is.null(side_by_side)) {
      side_by_side <- model_summary
    } else {
      side_by_side <- cbind(side_by_side, model_summary)
    }
  }
  return(side_by_side)
}

MODELLING Persuasion

# For indistinguishable Dyads
model_rows_fixed_persuasion_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_persuasion_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_persuasion_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:persuasion_self_cw',
    'is_B:persuasion_self_cw', 
    'is_A:persuasion_partner_cw', 
    'is_B:persuasion_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:persuasion_self_cb',
    'is_B:persuasion_self_cb', 
    'is_A:persuasion_partner_cb', 
    'is_B:persuasion_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_persuasion_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:persuasion_self_cw)',
  'sd(is_B:persuasion_self_cw)', 
  'sd(is_A:persuasion_partner_cw)', 
  'sd(is_B:persuasion_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ persuasion

range(df_double$pa_sub, na.rm = T) 

## [1]   0 720

hist(df_double$pa_sub, breaks = 100)

Modelling using the gaussian family fails. Due to the many zeros, transformations won’t help estimating the models. We employ the negative binomial family.

Indistinguishable

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day + 
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "persuasion_pa_sub_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_pa_sub_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_sub_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_pa_sub_ind)
print(loo_ind)

## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12082.9 177.5
## p_loo        28.2   2.4
## looic     24165.8 355.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.4]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_pa_sub_ind, ask = FALSE)

summarize_brms(
  persuasion_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()

	IRR	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	29.14*	20.18	41.57	1.001	903.07	1765.89
Within-Person Effects
persuasion_self_cw	1.01	0.92	1.11	1.000	9103.37	3111.55
persuasion_partner_cw	1.08	0.98	1.19	1.000	7783.86	2766.72
day	1.15	0.85	1.57	1.002	9017.24	3175.00
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
persuasion_self_cb	0.85	0.60	1.18	1.001	4385.12	2696.25
persuasion_partner_cb	1.11	0.79	1.59	1.000	3990.53	3056.64
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.72	0.52	0.98	1.01	1069.15	2159.42
sd(persuasion_self_cw)	0.25	0.12	0.40	1.00	2049.24	2247.12
sd(persuasion_partner_cw)	0.23	0.10	0.38	1.00	1924.16	2229.31
Additional Parameters
ar[1]	0.02	-0.94	0.94	1.00	6715.11	2833.16
nu	NA	NA	NA	NA	NA	NA
shape	0.14	0.13	0.14	1.00	8147.96	2612.93
sderr	0.05	0.00	0.13	1.00	2816.69	2191.39
sigma	NA	NA	NA	NA	NA	NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_pa_sub_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_pa_sub_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_sub_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_pa_sub_apim)

## Warning: Found 4 observations with a pareto_k > 0.7 in model 'persuasion_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12085.1 177.4
## p_loo        40.5   2.7
## looic     24170.2 354.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3738  99.9%   206     
##    (0.7, 1]   (bad)         4   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(persuasion_pa_sub_apim, ask = FALSE)

summarize_brms(
  persuasion_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()

	IRR	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	31.64*	21.06	47.91	1.001	1568.57	2509.33
is_B	20.52*	14.28	29.88	1.000	2136.80	2886.58
Within-Person Effects
is_A:persuasion_self_cw	0.98	0.86	1.14	1.002	7023.50	3392.76
is_B:persuasion_self_cw	1.05	0.91	1.21	1.001	6912.83	3178.03
is_A:persuasion_partner_cw	0.99	0.87	1.14	1.001	8085.20	2979.63
is_B:persuasion_partner_cw	1.17*	1.01	1.36	1.000	6556.89	3694.90
is_A:day	1.19	0.75	1.88	1.000	5424.27	2903.84
is_B:day	1.27	0.82	1.98	1.001	6204.54	3147.04
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:persuasion_self_cb	0.93	0.57	1.49	1.001	3673.30	3212.62
is_B:persuasion_self_cb	0.64	0.37	1.09	1.001	3619.34	3300.48
is_A:persuasion_partner_cb	1.09	0.65	1.87	1.000	3627.47	2952.14
is_B:persuasion_partner_cb	1.16	0.71	1.91	1.001	3721.53	3193.61
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.72	0.51	0.97	1.00	1950.14	2855.84
sd(is_B)	0.64	0.43	0.89	1.00	1661.64	2675.85
sd(is_A:persuasion_self_cw)	0.22	0.04	0.41	1.00	1392.37	1070.47
sd(is_B:persuasion_self_cw)	0.23	0.04	0.44	1.00	1545.44	1293.99
sd(is_A:persuasion_partner_cw)	0.18	0.01	0.40	1.00	1508.09	1470.95
sd(is_B:persuasion_partner_cw)	0.25	0.06	0.46	1.00	1662.67	1452.00
Additional Parameters
ar[1]	0.01	-0.94	0.96	1.00	5779.57	2085.73
nu	NA	NA	NA	NA	NA	NA
shape	0.14	0.13	0.14	1.00	6750.82	2677.34
sderr	0.04	0.00	0.13	1.00	2791.50	2417.53
sigma	NA	NA	NA	NA	NA	NA

Compare models

loo_compare(loo_ind, loo_apim)

##                        elpd_diff se_diff
## persuasion_pa_sub_ind   0.0       0.0   
## persuasion_pa_sub_apim -2.2       3.7

bayes_R2(persuasion_pa_sub_apim)

##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.2106389 0.0702662 0.1169052 0.4020931

bayes_R2(persuasion_pa_sub_ind)

##     Estimate  Est.Error       Q2.5     Q97.5
## R2 0.1703761 0.05028759 0.09973377 0.2920029

Device Based MVPA ~ persuasion

range(df_double$pa_obj, na.rm = T) 

## [1]   5.75 971.25

hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the model did not converge. Poisson also did not work. As we have no zeros in this distribution, we log transform.

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_pa_obj_log_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_pa_obj_log_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_obj_log_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_pa_obj_log_ind)
print(loo_ind)

## 
## Computed from 4000 by 3305 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2803.1  54.8
## p_loo        70.8   3.0
## looic      5606.2 109.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.5]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  persuasion_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()

## Warning in summarize_brms(persuasion_pa_obj_log_ind, model_rows_fixed = model_rows_fixed_persuasion_ind, :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	120.38*	108.19	134.27	1.011	458.70	1046.96
Within-Person Effects
persuasion_self_cw	1.01	0.99	1.03	1.001	6192.13	3094.03
persuasion_partner_cw	1.01	0.99	1.02	1.002	5939.04	3243.98
day	0.99	0.91	1.08	1.001	4877.50	3218.69
weartime_self_cw	1.00	1.00	1.00	1.002	3797.90	2901.27
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
persuasion_self_cb	0.87*	0.79	0.97	1.001	2456.65	2873.55
persuasion_partner_cb	1.05	0.95	1.17	1.001	2390.52	2787.94
weartime_self_cb	1.00	1.00	1.00	1.001	4350.93	3515.45
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.30	0.23	0.38	1.01	758.27	1077.55
sd(persuasion_self_cw)	0.06	0.03	0.09	1.00	2263.37	2401.22
sd(persuasion_partner_cw)	0.05	0.03	0.08	1.00	2075.45	1844.66
Additional Parameters
ar[1]	0.30	0.26	0.33	1.00	5720.76	2947.82
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.56	0.55	0.57	1.00	5237.61	3296.42

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_pa_obj_log_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_pa_obj_log_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_obj_log_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_pa_obj_log_apim)
print(loo_apim)

## 
## Computed from 4000 by 3305 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2721.8  56.3
## p_loo       129.9   5.7
## looic      5443.6 112.5
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.2]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  persuasion_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()

## Warning in summarize_brms(persuasion_pa_obj_log_apim, model_rows_fixed = model_rows_fixed_persuasion_apim,
## : Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	129.47*	108.37	153.25	1.005	675.45	1343.39
is_B	118.68*	101.47	140.04	1.003	717.81	1147.66
Within-Person Effects
is_A:persuasion_self_cw	0.97	0.95	1.00	1.001	6903.31	2865.10
is_B:persuasion_self_cw	1.05*	1.02	1.07	1.000	8393.27	3023.29
is_A:persuasion_partner_cw	1.00	0.97	1.03	1.000	10350.62	2548.48
is_B:persuasion_partner_cw	1.01	0.98	1.04	1.001	8481.91	2454.68
is_A:day	0.89*	0.79	0.99	1.000	5695.72	2639.66
is_B:day	1.06	0.95	1.19	1.004	6999.61	3316.51
is_A:weartime_self_cw	1.00	1.00	1.00	1.000	4098.01	2327.61
is_B:weartime_self_cw	1.00	1.00	1.00	1.001	3896.57	3064.60
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:persuasion_self_cb	0.89	0.77	1.04	1.000	2996.19	3204.21
is_B:persuasion_self_cb	0.85	0.71	1.00	1.000	2505.04	2772.26
is_A:persuasion_partner_cb	1.08	0.93	1.26	1.001	3333.22	2973.30
is_B:persuasion_partner_cb	1.15	0.99	1.34	1.001	2837.75	2698.05
is_A:weartime_self_cb	1.00	1.00	1.00	1.000	4489.33	3458.59
is_B:weartime_self_cb	1.00*	1.00	1.00	1.000	4627.69	3711.40
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.31	0.23	0.42	1.00	918.75	1454.42
sd(is_B)	0.39	0.29	0.52	1.00	736.28	1291.88
sd(is_A:persuasion_self_cw)	0.05	0.01	0.10	1.00	1396.31	1172.67
sd(is_B:persuasion_self_cw)	0.08	0.05	0.13	1.00	2783.30	3043.30
sd(is_A:persuasion_partner_cw)	0.07	0.03	0.13	1.00	1588.95	1861.31
sd(is_B:persuasion_partner_cw)	0.05	0.00	0.10	1.00	1064.17	842.20
Additional Parameters
ar[1]	0.22	0.18	0.25	1.00	6900.91	2930.32
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.54	0.53	0.55	1.00	7892.53	2411.23

Compare Models

loo_compare(loo_ind, loo_apim)

##                            elpd_diff se_diff
## persuasion_pa_obj_log_apim   0.0       0.0  
## persuasion_pa_obj_log_ind  -81.3      14.3

bayes_R2(persuasion_pa_obj_log_apim)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1890189 0.006237183 0.1766613 0.2013082

bayes_R2(persuasion_pa_obj_log_ind)

##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.1678555 0.006690044 0.1543916 0.180971

Affect ~ persuasion

range(df_double$aff, na.rm = T) 

## [1] 1 6

hist(df_double$aff, breaks = 15)

Indistinguishable

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_aff_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_aff_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_aff_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_aff_ind)
print(loo_ind)

## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4821.9  63.5
## p_loo        63.0   2.6
## looic      9643.9 126.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_aff_ind, ask = FALSE)

summarize_brms(
  persuasion_aff_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	4.90*	4.70	5.11	1.002	760.07	1412.86
Within-Person Effects
persuasion_self_cw	-0.01	-0.04	0.02	1.002	6135.29	2936.42
persuasion_partner_cw	0.03	0.00	0.06	1.001	7257.53	2884.58
day	-0.10	-0.25	0.04	1.000	6739.82	3092.48
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
persuasion_self_cb	0.02	-0.16	0.21	1.002	2929.97	2540.56
persuasion_partner_cb	0.17	-0.02	0.35	1.003	2577.32	2745.95
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.60	0.48	0.77	1.00	1006.48	1542.43
sd(persuasion_self_cw)	0.03	0.00	0.07	1.00	1663.48	1453.32
sd(persuasion_partner_cw)	0.07	0.01	0.12	1.00	929.44	1025.45
Additional Parameters
ar[1]	0.45	0.42	0.48	1.00	6242.15	2817.43
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.87	0.85	0.89	1.00	5910.04	2965.81

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_aff_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_aff_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_aff_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_aff_apim)
print(loo_apim)

## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4694.8  65.2
## p_loo       118.3   5.0
## looic      9389.6 130.5
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.4]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_aff_apim, ask = FALSE)

summarize_brms(
  persuasion_aff_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	4.88*	4.61	5.14	1.001	713.17	1297.44
is_B	4.89*	4.61	5.16	1.001	771.84	1579.80
Within-Person Effects
is_A:persuasion_self_cw	-0.02	-0.06	0.02	1.002	8436.28	2944.95
is_B:persuasion_self_cw	0.00	-0.04	0.04	1.000	9076.31	3247.40
is_A:persuasion_partner_cw	0.03	-0.01	0.07	1.003	7295.23	2741.16
is_B:persuasion_partner_cw	0.03	-0.01	0.07	1.002	7076.20	3073.96
is_A:day	-0.03	-0.21	0.16	1.004	6493.63	2849.65
is_B:day	-0.16	-0.34	0.02	1.003	8342.01	3166.59
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:persuasion_self_cb	0.05	-0.19	0.28	1.000	3092.29	3382.04
is_B:persuasion_self_cb	-0.12	-0.41	0.17	1.001	2558.03	2660.06
is_A:persuasion_partner_cb	0.07	-0.18	0.32	1.001	3230.63	3103.19
is_B:persuasion_partner_cb	0.01	-0.26	0.29	1.002	2399.79	2699.47
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.80	0.63	1.03	1.00	1146.74	1999.75
sd(is_B)	0.68	0.53	0.87	1.00	955.61	1858.38
sd(is_A:persuasion_self_cw)	0.05	0.00	0.11	1.01	1036.51	1603.46
sd(is_B:persuasion_self_cw)	0.08	0.01	0.15	1.00	1336.77	1537.21
sd(is_A:persuasion_partner_cw)	0.05	0.00	0.12	1.00	1388.60	1697.67
sd(is_B:persuasion_partner_cw)	0.08	0.01	0.15	1.00	1193.25	1533.01
Additional Parameters
ar[1]	0.34	0.30	0.37	1.00	5874.13	2813.52
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.83	0.82	0.85	1.00	7624.11	3099.01

Compare Models

loo_compare(loo_ind, loo_apim)

##                     elpd_diff se_diff
## persuasion_aff_apim    0.0       0.0 
## persuasion_aff_ind  -127.1      18.9

bayes_R2(persuasion_aff_apim)

##     Estimate  Est.Error    Q2.5     Q97.5
## R2 0.2341884 0.00451868 0.22503 0.2427378

bayes_R2(persuasion_aff_ind)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2247693 0.004741361 0.2155174 0.2337927

Anticipatory Affect ~ persuasion

range(df_double$anticipAff, na.rm = T) 

## [1] -5  5

hist(df_double$anticipAff, breaks = 10)

Indistinguishable

formula <- bf(
  anticipAff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_anticipAff_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_anticipAff_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_anticipAff_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_anticipAff_ind)
print(loo_ind)

## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7536.8  56.6
## p_loo        64.3   2.6
## looic     15073.6 113.2
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_anticipAff_ind, ask = FALSE)

summarize_brms(
  persuasion_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	1.47*	1.01	1.93	1.012	374.59	740.14
Within-Person Effects
persuasion_self_cw	-0.03	-0.08	0.03	1.000	4294.03	3090.91
persuasion_partner_cw	0.08*	0.02	0.14	1.000	4654.38	3167.04
day	-0.01	-0.32	0.30	1.001	4311.50	2389.89
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
persuasion_self_cb	-0.21	-0.57	0.15	1.002	1726.57	2121.35
persuasion_partner_cb	0.43*	0.05	0.80	1.002	1930.37	2210.76
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	1.27	0.99	1.60	1.01	583.64	1095.04
sd(persuasion_self_cw)	0.16	0.07	0.27	1.00	1388.45	1511.56
sd(persuasion_partner_cw)	0.04	0.00	0.11	1.00	1559.41	1549.23
Additional Parameters
ar[1]	0.42	0.39	0.45	1.00	3932.52	2916.81
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.80	1.76	1.84	1.00	5348.52	2742.81

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 5000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_anticipAff_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_anticipAff_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_anticipAff_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_anticipAff_apim)
print(loo_apim)

## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7439.2  56.7
## p_loo       114.1   4.2
## looic     14878.4 113.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_anticipAff_apim, ask = FALSE)

summarize_brms(
  persuasion_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	1.23*	0.60	1.86	1.004	607.10	1139.00
is_B	1.69*	1.15	2.27	1.002	817.88	1271.71
Within-Person Effects
is_A:persuasion_self_cw	-0.01	-0.09	0.07	1.003	8297.26	2853.71
is_B:persuasion_self_cw	-0.03	-0.11	0.05	1.003	8454.33	2857.25
is_A:persuasion_partner_cw	0.09*	0.01	0.17	1.001	8564.71	2885.68
is_B:persuasion_partner_cw	0.07	-0.01	0.15	1.000	8440.49	3254.58
is_A:day	0.06	-0.33	0.44	1.002	7256.56	2746.37
is_B:day	-0.08	-0.45	0.30	1.002	6818.28	3537.89
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:persuasion_self_cb	-0.22	-0.70	0.27	1.001	3493.49	3307.53
is_B:persuasion_self_cb	-0.59	-1.17	0.00	1.000	2676.13	2852.30
is_A:persuasion_partner_cb	0.43	-0.10	0.94	1.001	3815.72	3570.18
is_B:persuasion_partner_cb	0.50	-0.03	1.04	1.000	2770.49	2930.67
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	1.46	1.15	1.85	1.00	866.25	1926.05
sd(is_B)	1.49	1.16	1.91	1.00	1079.48	1532.73
sd(is_A:persuasion_self_cw)	0.13	0.01	0.29	1.00	1122.08	2193.39
sd(is_B:persuasion_self_cw)	0.20	0.07	0.34	1.00	1629.01	1530.04
sd(is_A:persuasion_partner_cw)	0.06	0.00	0.18	1.00	1975.80	2001.79
sd(is_B:persuasion_partner_cw)	0.06	0.00	0.17	1.00	2245.16	2141.55
Additional Parameters
ar[1]	0.33	0.30	0.37	1.00	7088.02	2876.23
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.74	1.70	1.78	1.00	8090.93	3063.62

Compare Models

loo_compare(loo_ind, loo_apim)

##                            elpd_diff se_diff
## persuasion_anticipAff_apim   0.0       0.0  
## persuasion_anticipAff_ind  -97.6      15.8

bayes_R2(persuasion_anticipAff_apim)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2193938 0.004568784 0.2103348 0.2284799

bayes_R2(persuasion_anticipAff_ind)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2105418 0.004928154 0.2010232 0.2202312

reactance ~ persuasion

range(df_double$reactance, na.rm = T) 

## [1] 0 5

hist(df_double$reactance, breaks = 6)

Indistinguishable

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_reactance_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_reactance_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_reactance_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_reactance_ind)
print(loo_ind)

## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1101.1 33.5
## p_loo        54.8  5.3
## looic      2202.1 66.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(persuasion_reactance_ind, ask = FALSE)

summarize_brms(
  persuasion_reactance_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.71*	0.43	0.98	1.002	1841.18	2594.40
Within-Person Effects
persuasion_self_cw	-0.01	-0.06	0.05	1.001	7277.61	2455.13
persuasion_partner_cw	-0.04	-0.10	0.02	1.002	9742.45	2931.45
day	-0.15	-0.45	0.13	1.001	7404.57	3125.69
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
persuasion_self_cb	0.02	-0.21	0.25	1.000	4189.26	2960.39
persuasion_partner_cb	0.07	-0.18	0.31	1.000	4112.40	2982.04
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.61	0.46	0.79	1.00	1971.40	2825.33
sd(persuasion_self_cw)	0.11	0.02	0.22	1.00	671.59	973.55
sd(persuasion_partner_cw)	0.05	0.00	0.12	1.00	1730.57	1918.83
Additional Parameters
ar[1]	0.07	-0.01	0.14	1.00	6223.85	2979.48
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.00	0.95	1.06	1.00	4549.32	2608.52

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_reactance_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(persuasion_reactance_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_reactance_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_reactance_apim)

## Warning: Found 3 observations with a pareto_k > 0.7 in model 'persuasion_reactance_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1082.0 33.1
## p_loo        90.5  8.1
## looic      2164.1 66.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     753   99.6%   83      
##    (0.7, 1]   (bad)        3    0.4%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(persuasion_reactance_apim, ask = FALSE)

summarize_brms(
  persuasion_reactance_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	0.40*	0.04	0.75	1.001	1679.69	2415.79
is_B	0.95*	0.58	1.33	1.001	2051.50	2436.77
Within-Person Effects
is_A:persuasion_self_cw	0.00	-0.08	0.07	1.000	4755.55	2853.10
is_B:persuasion_self_cw	-0.01	-0.08	0.06	1.002	6490.63	2909.35
is_A:persuasion_partner_cw	-0.02	-0.09	0.06	1.000	6363.01	3371.58
is_B:persuasion_partner_cw	-0.08	-0.16	0.01	1.001	5639.58	2773.14
is_A:day	0.18	-0.21	0.57	1.000	4082.36	3290.51
is_B:day	-0.37	-0.78	0.04	1.000	4530.83	3189.96
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:persuasion_self_cb	0.42*	0.09	0.76	1.001	2517.84	2853.79
is_B:persuasion_self_cb	-0.26	-0.67	0.12	1.002	2994.88	2948.22
is_A:persuasion_partner_cb	-0.53*	-0.94	-0.14	1.001	2581.08	2864.87
is_B:persuasion_partner_cb	0.35*	0.03	0.68	1.001	2781.56	2553.57
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.70	0.51	0.96	1.00	1408.17	2080.97
sd(is_B)	0.66	0.49	0.91	1.00	1783.02	2768.51
sd(is_A:persuasion_self_cw)	0.21	0.07	0.35	1.00	680.11	868.96
sd(is_B:persuasion_self_cw)	0.08	0.00	0.19	1.01	1029.42	1675.40
sd(is_A:persuasion_partner_cw)	0.07	0.00	0.18	1.00	1891.93	2446.46
sd(is_B:persuasion_partner_cw)	0.10	0.01	0.22	1.00	1083.22	945.22
Additional Parameters
ar[1]	-0.03	-0.12	0.05	1.00	3084.44	2569.49
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.95	0.89	1.00	1.00	3168.19	2657.30

Compare Models

loo_compare(loo_ind, loo_apim)

##                           elpd_diff se_diff
## persuasion_reactance_apim   0.0       0.0  
## persuasion_reactance_ind  -19.0       7.2

bayes_R2(persuasion_reactance_apim)

##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2386339 0.01018588 0.2179162 0.2579992

bayes_R2(persuasion_reactance_ind)

##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.2194376 0.0112548 0.1976425 0.2413073

Report models for Persuasion side by side

Indistinguishable Models

persuasion_all_ind_models <- report_side_by_side(
  persuasion_pa_sub_ind,
  persuasion_pa_obj_log_ind,
  persuasion_aff_ind,
  persuasion_anticipAff_ind,
  persuasion_reactance_ind,
  
  model_rows_random = model_rows_random_persuasion_ind,
  model_rows_fixed = model_rows_fixed_persuasion_ind
) 

## [1] "persuasion_pa_sub_ind"
## [1] "persuasion_pa_obj_log_ind"

## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.

## [1] "persuasion_aff_ind"
## [1] "persuasion_anticipAff_ind"
## [1] "persuasion_reactance_ind"

# pretty printing
persuasion_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()

	Subjective MVPA		Device-Based MVPA		Mood		Anticipatory Affect		Reactance
	IRR persuasion_pa_sub_ind	95% CI persuasion_pa_sub_ind	exp(Est.) persuasion_pa_obj_log_ind	95% CI persuasion_pa_obj_log_ind	b persuasion_aff_ind	95% CI persuasion_aff_ind	b persuasion_anticipAff_ind	95% CI persuasion_anticipAff_ind	b persuasion_reactance_ind	95% CI persuasion_reactance_ind
Intercept	29.14*	[20.18, 41.57]	120.38*	[108.19, 134.27]	4.90*	[ 4.70, 5.11]	1.47*	[ 1.01, 1.93]	0.71*	[ 0.43, 0.98]
Within-Person Effects
persuasion_self_cw	1.01	[ 0.92, 1.11]	1.01	[ 0.99, 1.03]	-0.01	[-0.04, 0.02]	-0.03	[-0.08, 0.03]	-0.01	[-0.06, 0.05]
persuasion_partner_cw	1.08	[ 0.98, 1.19]	1.01	[ 0.99, 1.02]	0.03	[ 0.00, 0.06]	0.08*	[ 0.02, 0.14]	-0.04	[-0.10, 0.02]
day	1.15	[ 0.85, 1.57]	0.99	[ 0.91, 1.08]	-0.10	[-0.25, 0.04]	-0.01	[-0.32, 0.30]	-0.15	[-0.45, 0.13]
weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
persuasion_self_cb	0.85	[ 0.60, 1.18]	0.87*	[ 0.79, 0.97]	0.02	[-0.16, 0.21]	-0.21	[-0.57, 0.15]	0.02	[-0.21, 0.25]
persuasion_partner_cb	1.11	[ 0.79, 1.59]	1.05	[ 0.95, 1.17]	0.17	[-0.02, 0.35]	0.43*	[ 0.05, 0.80]	0.07	[-0.18, 0.31]
weartime_self_cb	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.72	[ 0.52, 0.98]	0.30	[0.23, 0.38]	0.60	[0.48, 0.77]	1.27	[0.99, 1.60]	0.61	[ 0.46, 0.79]
sd(persuasion_self_cw)	0.25	[ 0.12, 0.40]	0.06	[0.03, 0.09]	0.03	[0.00, 0.07]	0.16	[0.07, 0.27]	0.11	[ 0.02, 0.22]
sd(persuasion_partner_cw)	0.23	[ 0.10, 0.38]	0.05	[0.03, 0.08]	0.07	[0.01, 0.12]	0.04	[0.00, 0.11]	0.05	[ 0.00, 0.12]
Additional Parameters
ar[1]	0.02	[-0.94, 0.94]	0.30	[0.26, 0.33]	0.45	[0.42, 0.48]	0.42	[0.39, 0.45]	0.07	[-0.01, 0.14]
nu	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
shape	0.14	[ 0.13, 0.14]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sderr	0.05	[ 0.00, 0.13]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sigma	NA	[ NA, NA]	0.56	[0.55, 0.57]	0.87	[0.85, 0.89]	1.80	[1.76, 1.84]	1.00	[ 0.95, 1.06]

# prepare export to excel
all_models <- list()
all_models$Persuasion_ind <- prepare_for_export(persuasion_all_ind_models)

APIMs

persuasion_all_apim_models <- report_side_by_side(
  persuasion_pa_sub_apim,
  persuasion_pa_obj_log_apim,
  persuasion_aff_apim,
  persuasion_anticipAff_apim,
  persuasion_reactance_apim,
  
  model_rows_random = model_rows_random_persuasion_apim,
  model_rows_fixed = model_rows_fixed_persuasion_apim
) 

## [1] "persuasion_pa_sub_apim"
## [1] "persuasion_pa_obj_log_apim"

## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.

## [1] "persuasion_aff_apim"
## [1] "persuasion_anticipAff_apim"
## [1] "persuasion_reactance_apim"

# pretty printing
persuasion_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()

	Subjective MVPA		Device-Based MVPA		Mood		Anticipatory Affect		Reactance
	IRR persuasion_pa_sub_apim	95% CI persuasion_pa_sub_apim	exp(Est.) persuasion_pa_obj_log_apim	95% CI persuasion_pa_obj_log_apim	b persuasion_aff_apim	95% CI persuasion_aff_apim	b persuasion_anticipAff_apim	95% CI persuasion_anticipAff_apim	b persuasion_reactance_apim	95% CI persuasion_reactance_apim
is_A	31.64*	[21.06, 47.91]	129.47*	[108.37, 153.25]	4.88*	[ 4.61, 5.14]	1.23*	[ 0.60, 1.86]	0.40*	[ 0.04, 0.75]
is_B	20.52*	[14.28, 29.88]	118.68*	[101.47, 140.04]	4.89*	[ 4.61, 5.16]	1.69*	[ 1.15, 2.27]	0.95*	[ 0.58, 1.33]
Within-Person Effects
is_A:persuasion_self_cw	0.98	[ 0.86, 1.14]	0.97	[ 0.95, 1.00]	-0.02	[-0.06, 0.02]	-0.01	[-0.09, 0.07]	0.00	[-0.08, 0.07]
is_B:persuasion_self_cw	1.05	[ 0.91, 1.21]	1.05*	[ 1.02, 1.07]	0.00	[-0.04, 0.04]	-0.03	[-0.11, 0.05]	-0.01	[-0.08, 0.06]
is_A:persuasion_partner_cw	0.99	[ 0.87, 1.14]	1.00	[ 0.97, 1.03]	0.03	[-0.01, 0.07]	0.09*	[ 0.01, 0.17]	-0.02	[-0.09, 0.06]
is_B:persuasion_partner_cw	1.17*	[ 1.01, 1.36]	1.01	[ 0.98, 1.04]	0.03	[-0.01, 0.07]	0.07	[-0.01, 0.15]	-0.08	[-0.16, 0.01]
is_A:day	1.19	[ 0.75, 1.88]	0.89*	[ 0.79, 0.99]	-0.03	[-0.21, 0.16]	0.06	[-0.33, 0.44]	0.18	[-0.21, 0.57]
is_B:day	1.27	[ 0.82, 1.98]	1.06	[ 0.95, 1.19]	-0.16	[-0.34, 0.02]	-0.08	[-0.45, 0.30]	-0.37	[-0.78, 0.04]
is_A:weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:persuasion_self_cb	0.93	[ 0.57, 1.49]	0.89	[ 0.77, 1.04]	0.05	[-0.19, 0.28]	-0.22	[-0.70, 0.27]	0.42*	[ 0.09, 0.76]
is_B:persuasion_self_cb	0.64	[ 0.37, 1.09]	0.85	[ 0.71, 1.00]	-0.12	[-0.41, 0.17]	-0.59	[-1.17, 0.00]	-0.26	[-0.67, 0.12]
is_A:persuasion_partner_cb	1.09	[ 0.65, 1.87]	1.08	[ 0.93, 1.26]	0.07	[-0.18, 0.32]	0.43	[-0.10, 0.94]	-0.53*	[-0.94, -0.14]
is_B:persuasion_partner_cb	1.16	[ 0.71, 1.91]	1.15	[ 0.99, 1.34]	0.01	[-0.26, 0.29]	0.50	[-0.03, 1.04]	0.35*	[ 0.03, 0.68]
is_A:weartime_self_cb	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	1.00*	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.72	[ 0.51, 0.97]	0.31	[0.23, 0.42]	0.80	[0.63, 1.03]	1.46	[1.15, 1.85]	0.70	[ 0.51, 0.96]
sd(is_B)	0.64	[ 0.43, 0.89]	0.39	[0.29, 0.52]	0.68	[0.53, 0.87]	1.49	[1.16, 1.91]	0.66	[ 0.49, 0.91]
sd(is_A:persuasion_self_cw)	0.22	[ 0.04, 0.41]	0.05	[0.01, 0.10]	0.05	[0.00, 0.11]	0.13	[0.01, 0.29]	0.21	[ 0.07, 0.35]
sd(is_B:persuasion_self_cw)	0.23	[ 0.04, 0.44]	0.08	[0.05, 0.13]	0.08	[0.01, 0.15]	0.20	[0.07, 0.34]	0.08	[ 0.00, 0.19]
sd(is_A:persuasion_partner_cw)	0.18	[ 0.01, 0.40]	0.07	[0.03, 0.13]	0.05	[0.00, 0.12]	0.06	[0.00, 0.18]	0.07	[ 0.00, 0.18]
sd(is_B:persuasion_partner_cw)	0.25	[ 0.06, 0.46]	0.05	[0.00, 0.10]	0.08	[0.01, 0.15]	0.06	[0.00, 0.17]	0.10	[ 0.01, 0.22]
Additional Parameters
ar[1]	0.01	[-0.94, 0.96]	0.22	[0.18, 0.25]	0.34	[0.30, 0.37]	0.33	[0.30, 0.37]	-0.03	[-0.12, 0.05]
nu	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
shape	0.14	[ 0.13, 0.14]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sderr	0.04	[ 0.00, 0.13]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sigma	NA	[ NA, NA]	0.54	[0.53, 0.55]	0.83	[0.82, 0.85]	1.74	[1.70, 1.78]	0.95	[ 0.89, 1.00]

# prepare export to excel
all_models$Persuasion_apim <- prepare_for_export(persuasion_all_apim_models)

Modelling Pressure

# For indistinguishable Dyads
model_rows_fixed_pressure_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'pressure_self_cb',
    'pressure_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_pressure_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_pressure_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:pressure_self_cw',
    'is_B:pressure_self_cw', 
    'is_A:pressure_partner_cw', 
    'is_B:pressure_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:pressure_self_cb',
    'is_B:pressure_self_cb', 
    'is_A:pressure_partner_cb', 
    'is_B:pressure_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_pressure_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:pressure_self_cw)',
  'sd(is_B:pressure_self_cw)', 
  'sd(is_A:pressure_partner_cw)', 
  'sd(is_B:pressure_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ pressure

Indistinguishable

formula <- bf(
  pa_sub ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day + 
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_sub_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_pa_sub_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_sub_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_pa_sub_ind)
print(loo_ind)

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12087.0 177.4
## p_loo        20.0   1.5
## looic     24173.9 354.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 2.5]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pressure_pa_sub_ind, ask = FALSE)

summarize_brms(
  pressure_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()

	IRR	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	24.44*	18.79	32.14	1.001	2684.70	3008.73
Within-Person Effects
pressure_self_cw	0.98	0.79	1.23	1.001	9687.42	2765.59
pressure_partner_cw	0.86	0.69	1.10	1.003	9495.13	2763.51
day	1.18	0.86	1.62	1.002	9666.04	2220.42
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
pressure_self_cb	1.09	0.66	1.81	1.000	5155.30	2888.05
pressure_partner_cb	0.89	0.55	1.47	1.000	5460.05	3143.76
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.63	0.47	0.84	1.00	1523.48	2255.08
sd(pressure_self_cw)	0.13	0.01	0.39	1.00	2670.03	2192.36
sd(pressure_partner_cw)	0.16	0.00	0.47	1.00	2200.86	1745.79
Additional Parameters
ar[1]	0.01	-0.93	0.94	1.00	7307.17	2816.12
nu	NA	NA	NA	NA	NA	NA
shape	0.13	0.12	0.14	1.00	9890.25	2099.07
sderr	0.04	0.00	0.13	1.00	2367.18	1497.39
sigma	NA	NA	NA	NA	NA	NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_sub_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_pa_sub_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_sub_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_pa_sub_apim)

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'pressure_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12086.2 177.4
## p_loo        29.0   2.0
## looic     24172.4 354.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3731  99.9%   476     
##    (0.7, 1]   (bad)         5   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pressure_pa_sub_apim, ask = FALSE)

summarize_brms(
  pressure_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()

	IRR	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	29.61*	20.61	41.74	1.003	1776.95	2047.69
is_B	19.98*	14.03	28.42	1.002	1672.04	2176.58
Within-Person Effects
is_A:pressure_self_cw	0.78	0.59	1.07	1.001	6725.39	3351.31
is_B:pressure_self_cw	1.26	0.91	1.88	1.001	3866.85	1938.48
is_A:pressure_partner_cw	1.02	0.76	1.45	1.001	5490.25	2781.88
is_B:pressure_partner_cw	0.63*	0.44	0.95	1.000	4440.11	2881.48
is_A:day	1.07	0.69	1.68	1.000	4162.17	2860.47
is_B:day	1.20	0.78	1.81	1.003	4004.45	2855.89
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:pressure_self_cb	1.08	0.49	2.36	1.000	1977.84	2075.41
is_B:pressure_self_cb	1.16	0.45	3.08	1.000	2129.31	2374.48
is_A:pressure_partner_cb	0.79	0.29	2.13	1.001	2090.44	2171.75
is_B:pressure_partner_cb	0.84	0.41	1.76	1.001	2526.13	2873.82
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.66	0.47	0.90	1.00	1972.83	2833.96
sd(is_B)	0.57	0.37	0.80	1.00	1870.38	2642.32
sd(is_A:pressure_self_cw)	0.21	0.01	0.70	1.00	2502.52	2167.92
sd(is_B:pressure_self_cw)	0.21	0.01	0.72	1.00	2768.34	2212.50
sd(is_A:pressure_partner_cw)	0.25	0.01	0.85	1.00	2796.74	2566.93
sd(is_B:pressure_partner_cw)	0.24	0.01	0.82	1.00	1578.33	1954.60
Additional Parameters
ar[1]	0.02	-0.94	0.94	1.00	4840.27	2464.21
nu	NA	NA	NA	NA	NA	NA
shape	0.13	0.13	0.14	1.00	4994.16	2779.02
sderr	0.04	0.00	0.12	1.00	2265.82	1795.20
sigma	NA	NA	NA	NA	NA	NA

Compare models

loo_compare(loo_ind, loo_apim)

##                      elpd_diff se_diff
## pressure_pa_sub_apim  0.0       0.0   
## pressure_pa_sub_ind  -0.7       2.6

bayes_R2(pressure_pa_sub_apim)

##     Estimate  Est.Error       Q2.5     Q97.5
## R2 0.1782304 0.09020382 0.09058281 0.4855923

bayes_R2(pressure_pa_sub_ind)

##     Estimate  Est.Error       Q2.5     Q97.5
## R2 0.1167618 0.03779368 0.06777029 0.1964773

Device Based MVPA ~ pressure

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_obj_log_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_pa_obj_log_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_obj_log_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_pa_obj_log_ind)
print(loo_ind)

## 
## Computed from 4000 by 3299 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2821.5  54.8
## p_loo        52.9   2.4
## looic      5643.0 109.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.4]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pressure_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  pressure_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()

## Warning in summarize_brms(pressure_pa_obj_log_ind, model_rows_fixed = model_rows_fixed_pressure_ind, :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	118.46*	106.30	131.66	1.009	535.38	1156.50
Within-Person Effects
pressure_self_cw	1.00	0.95	1.05	1.000	6377.09	3006.92
pressure_partner_cw	0.96	0.91	1.02	1.001	6449.17	2850.72
day	0.97	0.89	1.06	1.002	5721.89	2678.48
weartime_self_cw	1.00	1.00	1.00	1.001	3885.21	2144.43
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
pressure_self_cb	0.78*	0.67	0.93	1.002	2469.56	2730.16
pressure_partner_cb	1.18*	1.02	1.37	1.002	2278.16	2571.02
weartime_self_cb	1.00	1.00	1.00	1.001	3872.79	3061.44
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.29	0.22	0.37	1.01	814.16	1512.49
sd(pressure_self_cw)	0.05	0.00	0.14	1.00	1496.27	1841.42
sd(pressure_partner_cw)	0.03	0.00	0.10	1.00	2216.36	2241.51
Additional Parameters
ar[1]	0.28	0.25	0.32	1.00	5765.91	2848.84
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.56	0.55	0.58	1.00	6556.32	3074.86

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_obj_log_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_pa_obj_log_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_obj_log_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_pa_obj_log_apim)

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'pressure_pa_obj_log_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 3299 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2748.8  56.1
## p_loo       104.1   5.3
## looic      5497.7 112.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3294  99.8%   481     
##    (0.7, 1]   (bad)         5   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pressure_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  pressure_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()

## Warning in summarize_brms(pressure_pa_obj_log_apim, model_rows_fixed = model_rows_fixed_pressure_apim, :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	127.96*	106.77	150.24	1.017	531.49	857.08
is_B	123.39*	104.27	145.22	1.012	424.74	1029.31
Within-Person Effects
is_A:pressure_self_cw	1.01	0.94	1.09	1.000	7896.70	2787.06
is_B:pressure_self_cw	0.99	0.92	1.06	1.001	8265.48	2975.70
is_A:pressure_partner_cw	0.98	0.90	1.06	1.001	7960.95	3211.18
is_B:pressure_partner_cw	0.96	0.90	1.03	1.001	7788.09	2798.46
is_A:day	0.89*	0.79	0.99	0.999	6857.79	3004.78
is_B:day	1.01	0.91	1.13	1.002	7204.53	3132.25
is_A:weartime_self_cw	1.00	1.00	1.00	1.001	3671.36	2458.77
is_B:weartime_self_cw	1.00	1.00	1.00	1.001	3917.54	2748.69
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:pressure_self_cb	0.96	0.74	1.24	1.001	2660.19	2858.16
is_B:pressure_self_cb	0.57*	0.43	0.75	1.003	2662.74	3120.95
is_A:pressure_partner_cb	1.21	0.87	1.68	1.001	2634.87	3109.90
is_B:pressure_partner_cb	1.55*	1.25	1.92	1.002	2684.09	3001.05
is_A:weartime_self_cb	1.00	1.00	1.00	1.000	4690.54	3404.46
is_B:weartime_self_cb	1.00*	1.00	1.00	1.000	4172.07	3650.51
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.30	0.22	0.42	1.01	667.67	1367.32
sd(is_B)	0.39	0.29	0.52	1.01	519.76	1115.01
sd(is_A:pressure_self_cw)	0.05	0.00	0.15	1.00	2864.03	2463.77
sd(is_B:pressure_self_cw)	0.08	0.00	0.24	1.00	1698.22	1980.12
sd(is_A:pressure_partner_cw)	0.13	0.00	0.37	1.00	1112.61	1897.56
sd(is_B:pressure_partner_cw)	0.04	0.00	0.13	1.00	3045.19	1813.05
Additional Parameters
ar[1]	0.21	0.18	0.25	1.00	5907.45	2933.90
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.55	0.53	0.56	1.00	6026.96	2831.33

Compare Models

loo_compare(loo_ind, loo_apim)

##                          elpd_diff se_diff
## pressure_pa_obj_log_apim   0.0       0.0  
## pressure_pa_obj_log_ind  -72.7      14.3

bayes_R2(pressure_pa_obj_log_apim)

##    Estimate   Est.Error     Q2.5     Q97.5
## R2 0.183356 0.006415822 0.170636 0.1958512

bayes_R2(pressure_pa_obj_log_ind)

##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.1622604 0.006731342 0.1487327 0.175401

Affect ~ pressure

Indistinguishable

formula <- bf(
  aff ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_aff_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_aff_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_aff_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_aff_ind)
print(loo_ind)

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4821.6  63.5
## p_loo        58.1   3.3
## looic      9643.3 126.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.8]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pressure_aff_ind, ask = FALSE)

summarize_brms(
  pressure_aff_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	4.89*	4.67	5.10	1.001	745.62	1272.54
Within-Person Effects
pressure_self_cw	-0.02	-0.10	0.05	1.000	5382.92	2754.52
pressure_partner_cw	-0.01	-0.08	0.06	1.000	6190.02	2830.31
day	-0.12	-0.27	0.03	1.002	5330.47	2728.44
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
pressure_self_cb	0.12	-0.16	0.41	1.000	3197.40	2640.58
pressure_partner_cb	0.03	-0.24	0.32	1.001	3407.75	2855.72
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.60	0.47	0.78	1.01	788.16	1635.25
sd(pressure_self_cw)	0.10	0.00	0.26	1.00	832.42	1247.63
sd(pressure_partner_cw)	0.15	0.01	0.32	1.00	1122.34	1808.10
Additional Parameters
ar[1]	0.45	0.42	0.48	1.00	5392.27	2816.35
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.87	0.85	0.89	1.00	5374.82	2743.43

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_aff_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_aff_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_aff_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_aff_apim)

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'pressure_aff_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4694.2  65.3
## p_loo       101.0   4.4
## looic      9388.3 130.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 2.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  238     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pressure_aff_apim, ask = FALSE)

summarize_brms(
  pressure_aff_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	4.87*	4.58	5.13	1.002	618.64	1002.47
is_B	4.88*	4.60	5.17	1.003	605.34	891.66
Within-Person Effects
is_A:pressure_self_cw	-0.03	-0.12	0.08	1.000	7705.26	2999.35
is_B:pressure_self_cw	-0.05	-0.15	0.06	1.001	6987.03	3275.07
is_A:pressure_partner_cw	-0.03	-0.13	0.06	1.000	6252.92	2711.34
is_B:pressure_partner_cw	0.01	-0.09	0.11	1.002	6068.39	2747.39
is_A:day	-0.08	-0.25	0.10	1.004	6759.69	2583.51
is_B:day	-0.13	-0.31	0.05	1.002	6070.96	2812.16
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:pressure_self_cb	0.07	-0.34	0.47	1.000	2497.18	2873.84
is_B:pressure_self_cb	0.12	-0.33	0.56	1.001	2743.62	2687.95
is_A:pressure_partner_cb	-0.02	-0.51	0.48	1.000	2586.37	2966.02
is_B:pressure_partner_cb	-0.12	-0.48	0.24	1.000	2636.66	3049.75
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.79	0.62	1.01	1.00	903.34	2049.14
sd(is_B)	0.70	0.55	0.90	1.00	1001.21	1833.24
sd(is_A:pressure_self_cw)	0.10	0.00	0.31	1.00	1343.70	1872.04
sd(is_B:pressure_self_cw)	0.13	0.01	0.39	1.00	1767.18	2080.14
sd(is_A:pressure_partner_cw)	0.25	0.02	0.60	1.00	1117.76	1683.55
sd(is_B:pressure_partner_cw)	0.09	0.00	0.25	1.00	2168.54	1873.17
Additional Parameters
ar[1]	0.34	0.30	0.37	1.00	6535.64	2752.72
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.84	0.82	0.86	1.00	6234.22	2867.03

Compare Models

loo_compare(loo_ind, loo_apim)

##                   elpd_diff se_diff
## pressure_aff_apim    0.0       0.0 
## pressure_aff_ind  -127.5      18.6

bayes_R2(pressure_aff_apim)

##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.2349984 0.004573314 0.2260495 0.243759

bayes_R2(pressure_aff_ind)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2257293 0.004871448 0.2162224 0.2350853

Anticipatory Affect ~ pressure

Indistinguishable

formula <- bf(
  anticipAff ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_anticipAff_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_anticipAff_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_anticipAff_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_anticipAff_ind)
print(loo_ind)

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7525.4  56.6
## p_loo        60.7   2.8
## looic     15050.8 113.2
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pressure_anticipAff_ind, ask = FALSE)

summarize_brms(
  pressure_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	1.46*	1.02	1.91	1.005	603.42	1144.57
Within-Person Effects
pressure_self_cw	0.04	-0.11	0.20	1.003	6227.17	2541.02
pressure_partner_cw	0.01	-0.14	0.16	1.000	6249.71	2958.94
day	0.01	-0.29	0.30	1.001	6655.31	2963.37
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
pressure_self_cb	0.44	-0.10	0.96	1.001	2885.48	2768.50
pressure_partner_cb	-0.32	-0.85	0.23	1.000	3041.60	2824.84
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	1.27	0.99	1.61	1.01	836.58	1818.51
sd(pressure_self_cw)	0.31	0.12	0.56	1.00	2121.08	2136.77
sd(pressure_partner_cw)	0.30	0.10	0.53	1.00	2000.95	1516.93
Additional Parameters
ar[1]	0.42	0.39	0.45	1.00	5240.60	2632.30
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.80	1.76	1.84	1.00	6329.26	3065.50

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_anticipAff_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_anticipAff_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_anticipAff_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_anticipAff_apim)

## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pressure_anticipAff_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7424.0  56.4
## p_loo       106.8   4.5
## looic     14847.9 112.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3734  99.9%   390     
##    (0.7, 1]   (bad)         2   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pressure_anticipAff_apim, ask = FALSE)

summarize_brms(
  pressure_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	1.04*	0.43	1.62	1.005	821.47	1628.44
is_B	1.84*	1.26	2.44	1.006	689.95	1212.41
Within-Person Effects
is_A:pressure_self_cw	-0.04	-0.25	0.17	1.001	8818.85	2730.05
is_B:pressure_self_cw	0.16	-0.06	0.38	1.000	8690.77	2508.23
is_A:pressure_partner_cw	-0.10	-0.32	0.12	1.000	7399.40	3049.02
is_B:pressure_partner_cw	0.10	-0.10	0.31	1.003	8972.74	2544.64
is_A:day	0.22	-0.15	0.59	1.000	7710.82	3259.29
is_B:day	-0.17	-0.56	0.21	1.002	6770.78	2914.78
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:pressure_self_cb	0.11	-0.74	0.99	1.002	3186.28	2386.68
is_B:pressure_self_cb	-0.56	-1.49	0.40	1.001	3634.63	2900.67
is_A:pressure_partner_cb	0.03	-1.08	1.11	1.003	3222.52	2956.47
is_B:pressure_partner_cb	0.44	-0.32	1.20	1.000	3554.36	3010.47
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	1.47	1.15	1.90	1.00	946.11	1892.12
sd(is_B)	1.49	1.16	1.92	1.01	1025.47	1516.31
sd(is_A:pressure_self_cw)	0.41	0.08	0.78	1.00	1295.17	1095.07
sd(is_B:pressure_self_cw)	0.33	0.05	0.71	1.00	1594.05	1563.60
sd(is_A:pressure_partner_cw)	0.41	0.13	0.81	1.00	3051.13	2709.20
sd(is_B:pressure_partner_cw)	0.24	0.01	0.63	1.00	1319.27	1629.73
Additional Parameters
ar[1]	0.34	0.30	0.37	1.00	6037.29	2884.95
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.74	1.70	1.78	1.01	8425.43	2231.59

Compare Models

loo_compare(loo_ind, loo_apim)

##                          elpd_diff se_diff
## pressure_anticipAff_apim    0.0       0.0 
## pressure_anticipAff_ind  -101.4      16.0

bayes_R2(pressure_anticipAff_apim)

##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2202662 0.00455434 0.2110166 0.2288988

bayes_R2(pressure_anticipAff_ind)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2112457 0.004909094 0.2016515 0.2206921

reactance ~ pressure

Indistinguishable

formula <- bf(
  reactance ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  control = list(adapt_delta = 0.90),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_reactance_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_reactance_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_reactance_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_reactance_ind)

## Warning: Found 9 observations with a pareto_k > 0.7 in model 'pressure_reactance_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_ind)

## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1067.3 33.5
## p_loo        60.8  6.2
## looic      2134.7 67.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     747   98.8%   254     
##    (0.7, 1]   (bad)        8    1.1%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pressure_reactance_ind, ask = FALSE)

summarize_brms(
  pressure_reactance_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.56*	0.33	0.79	1.001	1102.00	2248.74
Within-Person Effects
pressure_self_cw	-0.02	-0.11	0.07	1.002	5289.22	3244.13
pressure_partner_cw	-0.06	-0.18	0.06	1.002	6071.65	3117.34
day	-0.13	-0.38	0.12	1.000	5474.68	3199.35
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
pressure_self_cb	0.10	-0.20	0.39	1.001	2425.72	2868.45
pressure_partner_cb	-0.05	-0.36	0.26	1.001	2329.16	2722.78
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.55	0.42	0.71	1.00	1575.28	2521.12
sd(pressure_self_cw)	0.44	0.27	0.66	1.00	1718.00	2543.78
sd(pressure_partner_cw)	0.17	0.01	0.51	1.00	964.20	1628.90
Additional Parameters
ar[1]	0.02	-0.06	0.10	1.00	4689.14	3011.29
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.95	0.90	1.00	1.00	4366.53	3027.93

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  control = list(adapt_delta = 0.90, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_reactance_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pressure_reactance_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_reactance_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_reactance_apim)

## Warning: Found 13 observations with a pareto_k > 0.7 in model 'pressure_reactance_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1055.8 33.8
## p_loo        79.3  7.5
## looic      2111.5 67.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     743   98.3%   155     
##    (0.7, 1]   (bad)       13    1.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pressure_reactance_apim, ask = FALSE)

summarize_brms(
  pressure_reactance_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	0.29*	0.00	0.57	1.002	2749.60	3190.70
is_B	0.77*	0.46	1.07	1.001	2297.08	2841.29
Within-Person Effects
is_A:pressure_self_cw	0.09	-0.04	0.22	1.001	8356.25	2795.65
is_B:pressure_self_cw	-0.12	-0.24	0.01	1.004	9228.00	3480.37
is_A:pressure_partner_cw	0.03	-0.16	0.23	1.003	7949.62	3269.05
is_B:pressure_partner_cw	-0.10	-0.25	0.04	1.003	10265.45	3046.31
is_A:day	0.12	-0.25	0.49	1.004	6959.77	3167.17
is_B:day	-0.25	-0.63	0.15	1.000	6353.28	3401.10
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:pressure_self_cb	0.51	-0.03	1.05	1.001	3906.50	3115.36
is_B:pressure_self_cb	-0.33	-0.96	0.28	1.001	3785.08	3021.86
is_A:pressure_partner_cb	-0.69	-1.40	0.01	1.000	3981.60	3159.78
is_B:pressure_partner_cb	0.33	-0.14	0.83	1.000	4266.32	3111.87
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.55	0.40	0.74	1.00	2046.58	3141.28
sd(is_B)	0.50	0.35	0.69	1.00	2747.64	3247.90
sd(is_A:pressure_self_cw)	0.42	0.23	0.70	1.00	2428.64	3141.31
sd(is_B:pressure_self_cw)	0.48	0.25	0.84	1.00	2588.51	3414.47
sd(is_A:pressure_partner_cw)	0.36	0.03	0.85	1.00	1720.32	1800.65
sd(is_B:pressure_partner_cw)	0.10	0.00	0.35	1.00	2722.46	2590.97
Additional Parameters
ar[1]	-0.02	-0.10	0.07	1.00	5172.64	3003.46
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.92	0.87	0.97	1.00	4939.79	3108.00

Compare Models

loo_compare(loo_ind, loo_apim)

##                         elpd_diff se_diff
## pressure_reactance_apim   0.0       0.0  
## pressure_reactance_ind  -11.6       6.1

bayes_R2(pressure_reactance_apim)

##    Estimate   Est.Error     Q2.5     Q97.5
## R2 0.243469 0.009813226 0.223816 0.2621822

bayes_R2(pressure_reactance_ind)

##     Estimate   Est.Error     Q2.5     Q97.5
## R2 0.2322911 0.009985695 0.212067 0.2511857

Report models for pressure side by side

Indistinguishable Models

pressure_all_ind_models <- report_side_by_side(
  pressure_pa_sub_ind,
  pressure_pa_obj_log_ind,
  pressure_aff_ind,
  pressure_anticipAff_ind,
  pressure_reactance_ind,
  
  model_rows_random = model_rows_random_pressure_ind,
  model_rows_fixed = model_rows_fixed_pressure_ind
) 

## [1] "pressure_pa_sub_ind"
## [1] "pressure_pa_obj_log_ind"

## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.

## [1] "pressure_aff_ind"
## [1] "pressure_anticipAff_ind"
## [1] "pressure_reactance_ind"

# pretty printing
pressure_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()

	Subjective MVPA		Device-Based MVPA		Mood		Anticipatory Affect		Reactance
	IRR pressure_pa_sub_ind	95% CI pressure_pa_sub_ind	exp(Est.) pressure_pa_obj_log_ind	95% CI pressure_pa_obj_log_ind	b pressure_aff_ind	95% CI pressure_aff_ind	b pressure_anticipAff_ind	95% CI pressure_anticipAff_ind	b pressure_reactance_ind	95% CI pressure_reactance_ind
Intercept	24.44*	[18.79, 32.14]	118.46*	[106.30, 131.66]	4.89*	[ 4.67, 5.10]	1.46*	[ 1.02, 1.91]	0.56*	[ 0.33, 0.79]
Within-Person Effects
pressure_self_cw	0.98	[ 0.79, 1.23]	1.00	[ 0.95, 1.05]	-0.02	[-0.10, 0.05]	0.04	[-0.11, 0.20]	-0.02	[-0.11, 0.07]
pressure_partner_cw	0.86	[ 0.69, 1.10]	0.96	[ 0.91, 1.02]	-0.01	[-0.08, 0.06]	0.01	[-0.14, 0.16]	-0.06	[-0.18, 0.06]
day	1.18	[ 0.86, 1.62]	0.97	[ 0.89, 1.06]	-0.12	[-0.27, 0.03]	0.01	[-0.29, 0.30]	-0.13	[-0.38, 0.12]
weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
barriers_self_cw	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
pressure_self_cb	1.09	[ 0.66, 1.81]	0.78*	[ 0.67, 0.93]	0.12	[-0.16, 0.41]	0.44	[-0.10, 0.96]	0.10	[-0.20, 0.39]
pressure_partner_cb	0.89	[ 0.55, 1.47]	1.18*	[ 1.02, 1.37]	0.03	[-0.24, 0.32]	-0.32	[-0.85, 0.23]	-0.05	[-0.36, 0.26]
weartime_self_cb	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
barriers_self_cb	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.63	[ 0.47, 0.84]	0.29	[0.22, 0.37]	0.60	[0.47, 0.78]	1.27	[0.99, 1.61]	0.55	[ 0.42, 0.71]
sd(pressure_self_cw)	0.13	[ 0.01, 0.39]	0.05	[0.00, 0.14]	0.10	[0.00, 0.26]	0.31	[0.12, 0.56]	0.44	[ 0.27, 0.66]
sd(pressure_partner_cw)	0.16	[ 0.00, 0.47]	0.03	[0.00, 0.10]	0.15	[0.01, 0.32]	0.30	[0.10, 0.53]	0.17	[ 0.01, 0.51]
Additional Parameters
ar[1]	0.01	[-0.93, 0.94]	0.28	[0.25, 0.32]	0.45	[0.42, 0.48]	0.42	[0.39, 0.45]	0.02	[-0.06, 0.10]
nu	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
shape	0.13	[ 0.12, 0.14]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sderr	0.04	[ 0.00, 0.13]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sigma	NA	[ NA, NA]	0.56	[0.55, 0.58]	0.87	[0.85, 0.89]	1.80	[1.76, 1.84]	0.95	[ 0.90, 1.00]

# prepare export to excel
all_models$Pressure_ind <- prepare_for_export(pressure_all_ind_models)

APIMs

pressure_all_apim_models <- report_side_by_side(
  pressure_pa_sub_apim,
  pressure_pa_obj_log_apim,
  pressure_aff_apim,
  pressure_anticipAff_apim,
  pressure_reactance_apim,
  
  model_rows_random = model_rows_random_pressure_apim,
  model_rows_fixed = model_rows_fixed_pressure_apim
) 

## [1] "pressure_pa_sub_apim"
## [1] "pressure_pa_obj_log_apim"

## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.

## [1] "pressure_aff_apim"

## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## [1] "pressure_anticipAff_apim"
## [1] "pressure_reactance_apim"

## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

# pretty printing
pressure_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()

	Subjective MVPA		Device-Based MVPA		Mood		Anticipatory Affect		Reactance
	IRR pressure_pa_sub_apim	95% CI pressure_pa_sub_apim	exp(Est.) pressure_pa_obj_log_apim	95% CI pressure_pa_obj_log_apim	b pressure_aff_apim	95% CI pressure_aff_apim	b pressure_anticipAff_apim	95% CI pressure_anticipAff_apim	b pressure_reactance_apim	95% CI pressure_reactance_apim
is_A	29.61*	[20.61, 41.74]	127.96*	[106.77, 150.24]	4.87*	[ 4.58, 5.13]	1.04*	[ 0.43, 1.62]	0.29*	[ 0.00, 0.57]
is_B	19.98*	[14.03, 28.42]	123.39*	[104.27, 145.22]	4.88*	[ 4.60, 5.17]	1.84*	[ 1.26, 2.44]	0.77*	[ 0.46, 1.07]
Within-Person Effects
is_A:pressure_self_cw	0.78	[ 0.59, 1.07]	1.01	[ 0.94, 1.09]	-0.03	[-0.12, 0.08]	-0.04	[-0.25, 0.17]	0.09	[-0.04, 0.22]
is_B:pressure_self_cw	1.26	[ 0.91, 1.88]	0.99	[ 0.92, 1.06]	-0.05	[-0.15, 0.06]	0.16	[-0.06, 0.38]	-0.12	[-0.24, 0.01]
is_A:pressure_partner_cw	1.02	[ 0.76, 1.45]	0.98	[ 0.90, 1.06]	-0.03	[-0.13, 0.06]	-0.10	[-0.32, 0.12]	0.03	[-0.16, 0.23]
is_B:pressure_partner_cw	0.63*	[ 0.44, 0.95]	0.96	[ 0.90, 1.03]	0.01	[-0.09, 0.11]	0.10	[-0.10, 0.31]	-0.10	[-0.25, 0.04]
is_A:day	1.07	[ 0.69, 1.68]	0.89*	[ 0.79, 0.99]	-0.08	[-0.25, 0.10]	0.22	[-0.15, 0.59]	0.12	[-0.25, 0.49]
is_B:day	1.20	[ 0.78, 1.81]	1.01	[ 0.91, 1.13]	-0.13	[-0.31, 0.05]	-0.17	[-0.56, 0.21]	-0.25	[-0.63, 0.15]
is_A:weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cw	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
is_A:pressure_self_cb	1.08	[ 0.49, 2.36]	0.96	[ 0.74, 1.24]	0.07	[-0.34, 0.47]	0.11	[-0.74, 0.99]	0.51	[-0.03, 1.05]
is_B:pressure_self_cb	1.16	[ 0.45, 3.08]	0.57*	[ 0.43, 0.75]	0.12	[-0.33, 0.56]	-0.56	[-1.49, 0.40]	-0.33	[-0.96, 0.28]
is_A:pressure_partner_cb	0.79	[ 0.29, 2.13]	1.21	[ 0.87, 1.68]	-0.02	[-0.51, 0.48]	0.03	[-1.08, 1.11]	-0.69	[-1.40, 0.01]
is_B:pressure_partner_cb	0.84	[ 0.41, 1.76]	1.55*	[ 1.25, 1.92]	-0.12	[-0.48, 0.24]	0.44	[-0.32, 1.20]	0.33	[-0.14, 0.83]
is_A:weartime_self_cb	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	1.00*	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
is_B:barriers_self_cb	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(is_A)	0.66	[ 0.47, 0.90]	0.30	[0.22, 0.42]	0.79	[0.62, 1.01]	1.47	[1.15, 1.90]	0.55	[ 0.40, 0.74]
sd(is_B)	0.57	[ 0.37, 0.80]	0.39	[0.29, 0.52]	0.70	[0.55, 0.90]	1.49	[1.16, 1.92]	0.50	[ 0.35, 0.69]
sd(is_A:pressure_self_cw)	0.21	[ 0.01, 0.70]	0.05	[0.00, 0.15]	0.10	[0.00, 0.31]	0.41	[0.08, 0.78]	0.42	[ 0.23, 0.70]
sd(is_B:pressure_self_cw)	0.21	[ 0.01, 0.72]	0.08	[0.00, 0.24]	0.13	[0.01, 0.39]	0.33	[0.05, 0.71]	0.48	[ 0.25, 0.84]
sd(is_A:pressure_partner_cw)	0.25	[ 0.01, 0.85]	0.13	[0.00, 0.37]	0.25	[0.02, 0.60]	0.41	[0.13, 0.81]	0.36	[ 0.03, 0.85]
sd(is_B:pressure_partner_cw)	0.24	[ 0.01, 0.82]	0.04	[0.00, 0.13]	0.09	[0.00, 0.25]	0.24	[0.01, 0.63]	0.10	[ 0.00, 0.35]
Additional Parameters
ar[1]	0.02	[-0.94, 0.94]	0.21	[0.18, 0.25]	0.34	[0.30, 0.37]	0.34	[0.30, 0.37]	-0.02	[-0.10, 0.07]
nu	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
shape	0.13	[ 0.13, 0.14]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sderr	0.04	[ 0.00, 0.12]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sigma	NA	[ NA, NA]	0.55	[0.53, 0.56]	0.84	[0.82, 0.86]	1.74	[1.70, 1.78]	0.92	[ 0.87, 0.97]

# prepare export to excel
all_models$Pressure_apim <- prepare_for_export(pressure_all_apim_models)

Modelling Pushing

# For indistinguishable Dyads
model_rows_fixed_pushing_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_pushing_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_pushing_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:pushing_self_cw',
    'is_B:pushing_self_cw', 
    'is_A:pushing_partner_cw', 
    'is_B:pushing_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:pushing_self_cb',
    'is_B:pushing_self_cb', 
    'is_A:pushing_partner_cb', 
    'is_B:pushing_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_pushing_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:pushing_self_cw)',
  'sd(is_B:pushing_self_cw)', 
  'sd(is_A:pushing_partner_cw)', 
  'sd(is_B:pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ pushing

Indistinguishable

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + 
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_sub_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_pa_sub_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_sub_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_pa_sub_ind)

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'pushing_pa_sub_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_ind)

## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -6316.5  67.5
## p_loo        33.5   2.8
## looic     12633.0 135.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1339  99.9%   242     
##    (0.7, 1]   (bad)         1   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pushing_pa_sub_ind, ask = FALSE)

summarize_brms(
  pushing_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()

	IRR	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	49.28*	37.64	64.23	1.001	1655.14	2441.46
Within-Person Effects
pushing_self_cw	1.00	0.90	1.11	1.000	5361.40	3118.25
pushing_partner_cw	1.02	0.92	1.13	1.000	5970.95	2971.70
day	0.86	0.64	1.15	1.001	6379.73	2842.66
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	1.02	0.95	1.08	1.000	7012.08	3036.44
Between-Person Effects
pushing_self_cb	0.59	0.34	1.03	1.000	3443.47	3192.15
pushing_partner_cb	1.50	0.91	2.49	1.000	5389.19	3197.07
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	1.03	0.90	1.17	1.001	3903.33	2745.51
Random Effects
sd(Intercept)	0.63	0.44	0.87	1.00	1301.04	2033.20
sd(pushing_self_cw)	0.31	0.13	0.51	1.00	1453.47	1985.91
sd(pushing_partner_cw)	0.28	0.11	0.48	1.00	1436.65	1793.03
Additional Parameters
ar[1]	0.04	-0.93	0.95	1.00	5529.16	2688.10
nu	NA	NA	NA	NA	NA	NA
shape	0.43	0.40	0.47	1.00	5112.46	3066.21
sderr	0.04	0.00	0.12	1.00	2385.84	2277.18
sigma	NA	NA	NA	NA	NA	NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_sub_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_pa_sub_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_sub_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_pa_sub_apim)

## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pushing_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -6323.5  67.5
## p_loo        43.1   2.8
## looic     12647.1 135.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1338  99.9%   292     
##    (0.7, 1]   (bad)         2   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pushing_pa_sub_apim, ask = FALSE)

summarize_brms(
  pushing_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()

## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

	IRR	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	64.43*	46.15	89.96	1.000	2298.52	2733.60
is_B	35.31*	25.01	49.86	1.000	2570.13	3174.92
Within-Person Effects
is_A:pushing_self_cw	1.05	0.92	1.21	1.001	7211.92	2859.27
is_B:pushing_self_cw	0.92	0.79	1.09	1.002	6161.35	3421.44
is_A:pushing_partner_cw	1.03	0.89	1.20	1.000	7086.77	3063.78
is_B:pushing_partner_cw	1.00	0.87	1.15	1.001	6245.30	3188.98
is_A:day	0.73	0.48	1.12	1.000	5634.86	3152.56
is_B:day	1.07	0.69	1.65	1.000	5717.33	3215.41
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	1.03	0.93	1.14	1.002	6587.34	2507.94
is_B:barriers_self_cw	0.99	0.91	1.08	1.000	7715.28	3127.70
Between-Person Effects
is_A:pushing_self_cb	0.47	0.19	1.15	1.000	3547.21	3588.53
is_B:pushing_self_cb	0.73	0.35	1.52	1.000	5000.77	3408.97
is_A:pushing_partner_cb	1.31	0.68	2.60	1.001	5162.93	3354.12
is_B:pushing_partner_cb	2.03	0.84	4.97	1.001	4122.35	3312.50
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	1.08	0.89	1.31	1.000	4317.48	3377.06
is_B:barriers_self_cb	1.04	0.82	1.30	1.001	3580.54	2922.44
Random Effects
sd(is_A)	0.63	0.39	0.97	1.00	1165.20	1894.07
sd(is_B)	0.47	0.29	0.68	1.00	1978.49	2586.86
sd(is_A:pushing_self_cw)	0.12	0.01	0.33	1.00	1640.62	2109.66
sd(is_B:pushing_self_cw)	0.46	0.14	0.85	1.00	1295.91	1465.64
sd(is_A:pushing_partner_cw)	0.43	0.11	0.84	1.00	1172.68	1373.02
sd(is_B:pushing_partner_cw)	0.11	0.01	0.30	1.00	2517.62	2903.67
Additional Parameters
ar[1]	0.04	-0.92	0.94	1.01	5063.36	2243.66
nu	NA	NA	NA	NA	NA	NA
shape	0.43	0.40	0.47	1.00	4753.53	2795.10
sderr	0.05	0.00	0.13	1.00	2691.83	2478.37
sigma	NA	NA	NA	NA	NA	NA

Compare models

loo_compare(loo_ind, loo_apim)

##                     elpd_diff se_diff
## pushing_pa_sub_ind   0.0       0.0   
## pushing_pa_sub_apim -7.0       3.6

bayes_R2(pushing_pa_sub_apim)

##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2654618 0.07201488 0.1666211 0.4672379

bayes_R2(pushing_pa_sub_ind)

##    Estimate  Est.Error      Q2.5     Q97.5
## R2 0.213283 0.04997866 0.1402257 0.3408443

Device Based MVPA ~ pushing

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_obj_log_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_pa_obj_log_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_obj_log_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_pa_obj_log_ind)
print(loo_ind)

## 
## Computed from 4000 by 1203 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -981.7 34.0
## p_loo        64.5  4.7
## looic      1963.4 68.0
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pushing_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  pushing_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()

## Warning in summarize_brms(pushing_pa_obj_log_ind, model_rows_fixed = model_rows_fixed_pushing_ind, :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	127.21*	111.99	144.82	1.000	1752.75	2358.30
Within-Person Effects
pushing_self_cw	1.00	0.97	1.04	1.001	7809.82	2778.98
pushing_partner_cw	0.99	0.96	1.03	1.000	7975.04	2950.37
day	1.00	0.88	1.14	1.000	7752.54	3200.40
weartime_self_cw	1.00	1.00	1.00	1.001	3820.73	2574.65
barriers_self_cw	0.98	0.96	1.01	1.001	7182.40	2846.49
Between-Person Effects
pushing_self_cb	0.82	0.64	1.07	1.002	4513.32	3468.94
pushing_partner_cb	0.86	0.66	1.12	1.001	4752.33	3505.51
weartime_self_cb	1.00	1.00	1.00	1.002	5376.81	3249.84
barriers_self_cb	0.99	0.93	1.05	1.002	4331.99	3156.88
Random Effects
sd(Intercept)	0.30	0.22	0.41	1.00	1217.79	2263.40
sd(pushing_self_cw)	0.10	0.04	0.17	1.00	1221.10	1899.51
sd(pushing_partner_cw)	0.04	0.00	0.11	1.00	1322.25	1982.92
Additional Parameters
ar[1]	0.27	0.21	0.33	1.00	5406.99	2539.64
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.53	0.51	0.55	1.00	6367.72	2955.61

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_obj_log_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_pa_obj_log_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_obj_log_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_pa_obj_log_apim)

## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pushing_pa_obj_log_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 1203 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -972.6 35.2
## p_loo       106.4  7.1
## looic      1945.3 70.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1201  99.8%   86      
##    (0.7, 1]   (bad)         2   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pushing_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  pushing_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()

## Warning in summarize_brms(pushing_pa_obj_log_apim, model_rows_fixed = model_rows_fixed_pushing_apim, :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	129.10*	106.65	155.78	1.003	831.15	1018.63
is_B	126.93*	104.96	153.25	1.004	910.22	1647.17
Within-Person Effects
is_A:pushing_self_cw	0.99	0.95	1.04	1.001	5241.24	2970.02
is_B:pushing_self_cw	1.02	0.96	1.08	1.000	5683.20	2956.87
is_A:pushing_partner_cw	1.03	0.97	1.09	1.000	4955.72	2932.87
is_B:pushing_partner_cw	0.98	0.93	1.02	1.000	6065.71	3072.27
is_A:day	1.08	0.91	1.28	1.000	4253.08	3276.88
is_B:day	0.94	0.78	1.12	1.000	3761.56	2977.15
is_A:weartime_self_cw	1.00	1.00	1.00	1.001	4018.46	2772.97
is_B:weartime_self_cw	1.00	1.00	1.00	1.002	4229.94	2751.28
is_A:barriers_self_cw	1.01	0.98	1.04	1.000	5530.29	3127.49
is_B:barriers_self_cw	0.98	0.95	1.00	1.000	5543.34	2635.67
Between-Person Effects
is_A:pushing_self_cb	0.81	0.56	1.19	1.000	3048.72	2968.96
is_B:pushing_self_cb	0.77	0.52	1.13	1.001	2956.11	2972.07
is_A:pushing_partner_cb	1.12	0.77	1.64	1.001	3227.02	3033.76
is_B:pushing_partner_cb	0.74	0.48	1.16	1.002	2533.24	2648.69
is_A:weartime_self_cb	1.00	1.00	1.00	1.000	3410.36	3401.54
is_B:weartime_self_cb	1.00	1.00	1.00	1.000	4222.92	3290.03
is_A:barriers_self_cb	1.07	0.99	1.17	1.001	2854.71	2630.65
is_B:barriers_self_cb	0.98	0.89	1.09	1.001	2917.54	3217.70
Random Effects
sd(is_A)	0.34	0.24	0.47	1.00	1058.13	2038.29
sd(is_B)	0.34	0.23	0.48	1.00	857.30	1343.85
sd(is_A:pushing_self_cw)	0.09	0.01	0.18	1.01	911.24	911.67
sd(is_B:pushing_self_cw)	0.11	0.01	0.20	1.00	922.39	638.42
sd(is_A:pushing_partner_cw)	0.11	0.02	0.20	1.00	1053.84	885.89
sd(is_B:pushing_partner_cw)	0.03	0.00	0.09	1.00	2359.01	2068.39
Additional Parameters
ar[1]	0.21	0.14	0.28	1.00	3704.50	3396.22
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.52	0.50	0.54	1.00	3544.78	2792.99

Compare Models

loo_compare(loo_ind, loo_apim)

##                         elpd_diff se_diff
## pushing_pa_obj_log_apim  0.0       0.0   
## pushing_pa_obj_log_ind  -9.1       7.6

bayes_R2(pushing_pa_obj_log_apim)

##    Estimate  Est.Error     Q2.5     Q97.5
## R2 0.210763 0.01165598 0.187023 0.2326488

bayes_R2(pushing_pa_obj_log_ind)

##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.1859197  0.012287 0.1615149 0.2094954

Affect ~ pushing

Indistinguishable

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_aff_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_aff_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_aff_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_aff_ind)
print(loo_ind)

## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1769.7 41.4
## p_loo        55.0  3.4
## looic      3539.5 82.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pushing_aff_ind, ask = FALSE)

summarize_brms(
  pushing_aff_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	5.00*	4.77	5.24	1.008	705.32	1279.78
Within-Person Effects
pushing_self_cw	-0.03	-0.08	0.03	1.004	5478.74	3016.82
pushing_partner_cw	0.03	-0.02	0.09	1.001	5545.96	2960.22
day	-0.04	-0.24	0.15	1.001	4709.38	2754.74
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	0.00	-0.04	0.03	1.001	5371.78	3127.61
Between-Person Effects
pushing_self_cb	-0.07	-0.49	0.33	1.000	2734.11	3086.61
pushing_partner_cb	0.22	-0.16	0.60	1.000	3505.84	2982.04
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	-0.08	-0.18	0.02	1.001	3245.62	2991.88
Random Effects
sd(Intercept)	0.62	0.47	0.81	1.00	971.33	1624.25
sd(pushing_self_cw)	0.07	0.01	0.15	1.00	1114.99	1379.33
sd(pushing_partner_cw)	0.12	0.03	0.22	1.00	1204.91	725.56
Additional Parameters
ar[1]	0.27	0.22	0.33	1.00	4812.62	3064.25
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.89	0.86	0.92	1.00	5261.68	2709.15

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_aff_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_aff_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_aff_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_aff_apim)
print(loo_apim)

## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1740.0 41.3
## p_loo        98.2  5.7
## looic      3480.0 82.6
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pushing_aff_apim, ask = FALSE)

summarize_brms(
  pushing_aff_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	5.00*	4.71	5.29	1.001	1540.42	2422.92
is_B	4.90*	4.59	5.21	1.001	1294.81	2299.76
Within-Person Effects
is_A:pushing_self_cw	-0.01	-0.08	0.07	1.001	7720.22	2873.36
is_B:pushing_self_cw	-0.06	-0.14	0.03	1.000	8447.37	3343.11
is_A:pushing_partner_cw	0.00	-0.08	0.09	1.000	7345.86	2894.63
is_B:pushing_partner_cw	0.05	-0.02	0.13	1.002	8307.39	3042.29
is_A:day	-0.06	-0.33	0.21	1.000	6880.69	3113.49
is_B:day	-0.01	-0.26	0.25	1.001	7737.84	3345.38
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	-0.02	-0.07	0.04	1.003	8714.70	2540.72
is_B:barriers_self_cw	0.02	-0.03	0.06	1.001	7904.22	2670.31
Between-Person Effects
is_A:pushing_self_cb	0.42	-0.22	1.08	1.002	3505.56	3356.73
is_B:pushing_self_cb	-0.28	-0.84	0.26	1.000	4412.31	3472.92
is_A:pushing_partner_cb	-0.07	-0.61	0.49	1.000	4266.48	3331.04
is_B:pushing_partner_cb	0.41	-0.23	1.06	1.000	4079.82	3420.17
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	0.01	-0.14	0.16	1.000	3759.93	3437.66
is_B:barriers_self_cb	-0.07	-0.25	0.11	1.000	2812.54	2230.19
Random Effects
sd(is_A)	0.79	0.61	1.02	1.00	1664.25	2441.20
sd(is_B)	0.68	0.51	0.89	1.00	1527.28	2451.03
sd(is_A:pushing_self_cw)	0.07	0.00	0.18	1.00	2328.67	2045.82
sd(is_B:pushing_self_cw)	0.12	0.01	0.27	1.00	1568.23	1969.97
sd(is_A:pushing_partner_cw)	0.17	0.04	0.32	1.00	1934.22	1642.01
sd(is_B:pushing_partner_cw)	0.10	0.01	0.24	1.00	1728.04	2285.82
Additional Parameters
ar[1]	0.18	0.12	0.24	1.00	4212.94	2703.09
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	0.85	0.82	0.89	1.00	5806.11	2871.44

Compare Models

loo_compare(loo_ind, loo_apim)

##                  elpd_diff se_diff
## pushing_aff_apim   0.0       0.0  
## pushing_aff_ind  -29.7      10.3

bayes_R2(pushing_aff_apim)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2077469 0.007994099 0.1918707 0.2232979

bayes_R2(pushing_aff_ind)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1898545 0.008693353 0.1725112 0.2062479

Anticipatory Affect ~ pushing

Indistinguishable

formula <- bf(
  anticipAff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + 
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_anticipAff_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_anticipAff_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_anticipAff_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_anticipAff_ind)
print(loo_ind)

## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2770.8 30.9
## p_loo        51.8  2.6
## looic      5541.5 61.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.6]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

plot(pushing_anticipAff_ind, ask = FALSE)

summarize_brms(
  pushing_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	1.91*	1.44	2.37	1.007	804.61	1509.41
Within-Person Effects
pushing_self_cw	0.01	-0.11	0.12	1.001	5355.07	2643.69
pushing_partner_cw	-0.07	-0.19	0.05	1.003	4170.72	3032.32
day	-0.25	-0.66	0.19	1.000	4744.55	3205.48
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	-0.01	-0.09	0.06	1.000	5164.31	2685.59
Between-Person Effects
pushing_self_cb	-0.34	-1.22	0.55	1.001	2616.54	2962.99
pushing_partner_cb	-0.05	-0.91	0.82	1.000	3466.69	2656.87
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	0.17	-0.06	0.39	1.002	2864.71	2851.86
Random Effects
sd(Intercept)	1.22	0.92	1.60	1.00	1047.94	2128.53
sd(pushing_self_cw)	0.15	0.01	0.34	1.00	1188.27	1230.81
sd(pushing_partner_cw)	0.17	0.01	0.39	1.00	1192.06	1434.02
Additional Parameters
ar[1]	0.34	0.28	0.39	1.00	4137.46	3139.97
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.88	1.80	1.95	1.00	4819.75	2935.38

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_anticipAff_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_anticipAff_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_anticipAff_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_anticipAff_apim)

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'pushing_anticipAff_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2738.5 31.8
## p_loo        89.1  4.4
## looic      5477.0 63.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1339  99.9%   737     
##    (0.7, 1]   (bad)         1   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pushing_anticipAff_apim, ask = FALSE)

summarize_brms(
  pushing_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	1.60*	0.96	2.23	1.000	2336.15	2468.72
is_B	2.02*	1.41	2.65	1.002	2400.99	2975.07
Within-Person Effects
is_A:pushing_self_cw	0.04	-0.11	0.19	1.000	8718.69	3060.17
is_B:pushing_self_cw	-0.07	-0.25	0.10	1.001	8885.25	3221.91
is_A:pushing_partner_cw	0.00	-0.18	0.17	1.001	9114.76	2944.38
is_B:pushing_partner_cw	-0.09	-0.25	0.06	1.001	7658.91	3226.08
is_A:day	0.06	-0.51	0.64	1.001	7604.15	3302.30
is_B:day	-0.55	-1.13	0.05	1.000	9013.42	2681.08
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	0.01	-0.10	0.12	1.000	10377.14	3052.80
is_B:barriers_self_cw	-0.03	-0.13	0.06	1.001	8435.68	2823.84
Between-Person Effects
is_A:pushing_self_cb	1.01	-0.47	2.54	1.000	4815.92	3066.66
is_B:pushing_self_cb	-0.80	-1.96	0.41	1.002	5833.57	3178.21
is_A:pushing_partner_cb	-1.21	-2.51	0.05	1.000	5793.34	3572.47
is_B:pushing_partner_cb	1.16	-0.36	2.58	1.001	5095.68	3114.79
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	0.21	-0.11	0.54	1.000	5009.92	3312.75
is_B:barriers_self_cb	0.04	-0.35	0.42	1.001	4198.81	3167.66
Random Effects
sd(is_A)	1.43	1.08	1.88	1.00	1466.89	2692.09
sd(is_B)	1.52	1.15	1.98	1.00	2066.01	2792.50
sd(is_A:pushing_self_cw)	0.15	0.01	0.37	1.00	2017.59	2134.81
sd(is_B:pushing_self_cw)	0.19	0.01	0.49	1.00	1942.39	2438.13
sd(is_A:pushing_partner_cw)	0.28	0.02	0.60	1.00	1525.13	1808.11
sd(is_B:pushing_partner_cw)	0.11	0.00	0.31	1.00	2727.86	2989.60
Additional Parameters
ar[1]	0.24	0.18	0.30	1.00	5413.34	2912.22
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.80	1.73	1.88	1.00	7463.90	2647.65

Compare Models

loo_compare(loo_ind, loo_apim)

##                         elpd_diff se_diff
## pushing_anticipAff_apim   0.0       0.0  
## pushing_anticipAff_ind  -32.3       9.7

bayes_R2(pushing_anticipAff_apim)

##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.2212212 0.008630534 0.2035636 0.237974

bayes_R2(pushing_anticipAff_ind)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2016732 0.009264627 0.1836161 0.2195859

reactance ~ pushing

Indistinguishable

formula <- bf(
  reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_reactance_ind")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_reactance_ind, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_reactance_ind)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_reactance_ind)

## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pushing_reactance_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_ind)

## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -597.3 24.5
## p_loo        38.4  4.8
## looic      1194.6 49.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     401   99.5%   301     
##    (0.7, 1]   (bad)        2    0.5%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pushing_reactance_ind, ask = FALSE)

summarize_brms(
  pushing_reactance_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.48*	0.25	0.71	1.002	6101.16	3321.85
Within-Person Effects
pushing_self_cw	-0.02	-0.10	0.07	1.002	7139.75	3154.64
pushing_partner_cw	0.02	-0.07	0.11	1.001	8174.34	3071.49
day	0.27	-0.12	0.66	1.003	8161.51	2643.71
weartime_self_cw	NA	NA	NA	NA	NA	NA
barriers_self_cw	0.05	-0.03	0.12	1.003	8452.29	3346.46
Between-Person Effects
pushing_self_cb	-0.71*	-1.22	-0.20	1.000	5670.61	3274.64
pushing_partner_cb	0.27	-0.31	0.87	1.000	5673.75	3392.01
weartime_self_cb	NA	NA	NA	NA	NA	NA
barriers_self_cb	0.01	-0.17	0.18	1.002	4772.32	3303.19
Random Effects
sd(Intercept)	0.17	0.01	0.38	1.00	1653.80	2222.19
sd(pushing_self_cw)	0.23	0.11	0.37	1.00	1449.65	1987.04
sd(pushing_partner_cw)	0.08	0.00	0.23	1.00	1770.45	1884.97
Additional Parameters
ar[1]	0.16	0.03	0.28	1.00	5132.70	3219.90
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.01	0.94	1.09	1.00	5895.25	3394.64

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_reactance_apim")
)

## Warning: Rows containing NAs were excluded from the model.

pp_check(pushing_reactance_apim, type='hist')

## Using 10 posterior draws for ppc type 'hist' by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_reactance_apim)

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_reactance_apim)

## Warning: Found 4 observations with a pareto_k > 0.7 in model 'pushing_reactance_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.

print(loo_apim)

## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -603.3 24.4
## p_loo        58.7  6.8
## looic      1206.6 48.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     399   99.0%   118     
##    (0.7, 1]   (bad)        4    1.0%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(pushing_reactance_apim, ask = FALSE)

summarize_brms(
  pushing_reactance_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()

	b	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
is_A	0.30	-0.01	0.60	1.001	3000.47	2975.13
is_B	0.74*	0.41	1.08	1.004	2659.02	2802.83
Within-Person Effects
is_A:pushing_self_cw	0.00	-0.12	0.11	1.002	4557.74	3278.77
is_B:pushing_self_cw	-0.05	-0.18	0.08	1.000	4887.79	3100.72
is_A:pushing_partner_cw	0.00	-0.13	0.13	1.001	5093.91	3230.72
is_B:pushing_partner_cw	0.05	-0.09	0.18	1.001	4199.18	2933.38
is_A:day	0.30	-0.24	0.85	1.002	3657.22	3219.69
is_B:day	0.03	-0.58	0.62	1.002	2777.06	3062.17
is_A:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	0.08	-0.03	0.19	1.000	5788.53	2996.24
is_B:barriers_self_cw	0.01	-0.10	0.11	1.001	4903.70	2641.80
Between-Person Effects
is_A:pushing_self_cb	-0.40	-1.15	0.37	1.000	4213.05	3065.63
is_B:pushing_self_cb	-0.76	-1.52	0.04	1.001	2795.76	2872.61
is_A:pushing_partner_cb	0.10	-0.63	0.86	1.001	4719.23	2843.64
is_B:pushing_partner_cb	0.51	-0.54	1.59	1.000	2620.35	2881.72
is_A:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	0.10	-0.11	0.32	1.000	3490.33	2884.67
is_B:barriers_self_cb	-0.16	-0.46	0.14	1.000	4182.24	3300.63
Random Effects
sd(is_A)	0.16	0.01	0.42	1.00	1611.07	2199.42
sd(is_B)	0.26	0.01	0.64	1.01	980.82	1268.92
sd(is_A:pushing_self_cw)	0.29	0.11	0.48	1.00	1189.28	1251.93
sd(is_B:pushing_self_cw)	0.23	0.04	0.43	1.01	895.24	1039.10
sd(is_A:pushing_partner_cw)	0.13	0.00	0.39	1.00	917.83	1897.62
sd(is_B:pushing_partner_cw)	0.12	0.01	0.35	1.01	1445.05	1996.32
Additional Parameters
ar[1]	0.13	0.00	0.27	1.00	2502.46	2816.40
nu	NA	NA	NA	NA	NA	NA
shape	NA	NA	NA	NA	NA	NA
sderr	NA	NA	NA	NA	NA	NA
sigma	1.00	0.92	1.08	1.00	2258.23	2754.95

Compare Models

loo_compare(loo_ind, loo_apim)

##                        elpd_diff se_diff
## pushing_reactance_ind   0.0       0.0   
## pushing_reactance_apim -6.0       3.5

bayes_R2(pushing_reactance_apim)

##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.1868663 0.01869625 0.1482543 0.2219012

bayes_R2(pushing_reactance_ind)

##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.1611657 0.0191676 0.1218981 0.1972937

Report models for pushing side by side

Indistinguishable Models

pushing_all_ind_models <- report_side_by_side(
  pushing_pa_sub_ind,
  pushing_pa_obj_log_ind,
  pushing_aff_ind,
  pushing_anticipAff_ind,
  pushing_reactance_ind,
  
  model_rows_random = model_rows_random_pushing_ind,
  model_rows_fixed = model_rows_fixed_pushing_ind
) 

## [1] "pushing_pa_sub_ind"
## [1] "pushing_pa_obj_log_ind"

## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.

## [1] "pushing_aff_ind"
## [1] "pushing_anticipAff_ind"
## [1] "pushing_reactance_ind"

# pretty printing
pushing_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()

	Subjective MVPA		Device-Based MVPA		Mood		Anticipatory Affect		Reactance
	IRR pushing_pa_sub_ind	95% CI pushing_pa_sub_ind	exp(Est.) pushing_pa_obj_log_ind	95% CI pushing_pa_obj_log_ind	b pushing_aff_ind	95% CI pushing_aff_ind	b pushing_anticipAff_ind	95% CI pushing_anticipAff_ind	b pushing_reactance_ind	95% CI pushing_reactance_ind
Intercept	49.28*	[37.64, 64.23]	127.21*	[111.99, 144.82]	5.00*	[ 4.77, 5.24]	1.91*	[ 1.44, 2.37]	0.48*	[ 0.25, 0.71]
Within-Person Effects
pushing_self_cw	1.00	[ 0.90, 1.11]	1.00	[ 0.97, 1.04]	-0.03	[-0.08, 0.03]	0.01	[-0.11, 0.12]	-0.02	[-0.10, 0.07]
pushing_partner_cw	1.02	[ 0.92, 1.13]	0.99	[ 0.96, 1.03]	0.03	[-0.02, 0.09]	-0.07	[-0.19, 0.05]	0.02	[-0.07, 0.11]
day	0.86	[ 0.64, 1.15]	1.00	[ 0.88, 1.14]	-0.04	[-0.24, 0.15]	-0.25	[-0.66, 0.19]	0.27	[-0.12, 0.66]
weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
barriers_self_cw	1.02	[ 0.95, 1.08]	0.98	[ 0.96, 1.01]	0.00	[-0.04, 0.03]	-0.01	[-0.09, 0.06]	0.05	[-0.03, 0.12]
Between-Person Effects
pushing_self_cb	0.59	[ 0.34, 1.03]	0.82	[ 0.64, 1.07]	-0.07	[-0.49, 0.33]	-0.34	[-1.22, 0.55]	-0.71*	[-1.22, -0.20]
pushing_partner_cb	1.50	[ 0.91, 2.49]	0.86	[ 0.66, 1.12]	0.22	[-0.16, 0.60]	-0.05	[-0.91, 0.82]	0.27	[-0.31, 0.87]
weartime_self_cb	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
barriers_self_cb	1.03	[ 0.90, 1.17]	0.99	[ 0.93, 1.05]	-0.08	[-0.18, 0.02]	0.17	[-0.06, 0.39]	0.01	[-0.17, 0.18]
Random Effects
sd(Intercept)	0.63	[ 0.44, 0.87]	0.30	[0.22, 0.41]	0.62	[0.47, 0.81]	1.22	[0.92, 1.60]	0.17	[0.01, 0.38]
sd(pushing_self_cw)	0.31	[ 0.13, 0.51]	0.10	[0.04, 0.17]	0.07	[0.01, 0.15]	0.15	[0.01, 0.34]	0.23	[0.11, 0.37]
sd(pushing_partner_cw)	0.28	[ 0.11, 0.48]	0.04	[0.00, 0.11]	0.12	[0.03, 0.22]	0.17	[0.01, 0.39]	0.08	[0.00, 0.23]
Additional Parameters
ar[1]	0.04	[-0.93, 0.95]	0.27	[0.21, 0.33]	0.27	[0.22, 0.33]	0.34	[0.28, 0.39]	0.16	[0.03, 0.28]
nu	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
shape	0.43	[ 0.40, 0.47]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sderr	0.04	[ 0.00, 0.12]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sigma	NA	[ NA, NA]	0.53	[0.51, 0.55]	0.89	[0.86, 0.92]	1.88	[1.80, 1.95]	1.01	[0.94, 1.09]

# prepare export to excel
all_models$Pushing_ind <- prepare_for_export(pushing_all_ind_models)

APIMs

pushing_all_apim_models <- report_side_by_side(
  pushing_pa_sub_apim,
  pushing_pa_obj_log_apim,
  pushing_aff_apim,
  pushing_anticipAff_apim,
  pushing_reactance_apim,
  
  model_rows_random = model_rows_random_pushing_apim,
  model_rows_fixed = model_rows_fixed_pushing_apim
) 

## [1] "pushing_pa_sub_apim"

## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## [1] "pushing_pa_obj_log_apim"

## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.

## [1] "pushing_aff_apim"
## [1] "pushing_anticipAff_apim"
## [1] "pushing_reactance_apim"

# pretty printing
pushing_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()

	Subjective MVPA		Device-Based MVPA		Mood		Anticipatory Affect		Reactance
	IRR pushing_pa_sub_apim	95% CI pushing_pa_sub_apim	exp(Est.) pushing_pa_obj_log_apim	95% CI pushing_pa_obj_log_apim	b pushing_aff_apim	95% CI pushing_aff_apim	b pushing_anticipAff_apim	95% CI pushing_anticipAff_apim	b pushing_reactance_apim	95% CI pushing_reactance_apim
is_A	64.43*	[46.15, 89.96]	129.10*	[106.65, 155.78]	5.00*	[ 4.71, 5.29]	1.60*	[ 0.96, 2.23]	0.30	[-0.01, 0.60]
is_B	35.31*	[25.01, 49.86]	126.93*	[104.96, 153.25]	4.90*	[ 4.59, 5.21]	2.02*	[ 1.41, 2.65]	0.74*	[ 0.41, 1.08]
Within-Person Effects
is_A:pushing_self_cw	1.05	[ 0.92, 1.21]	0.99	[ 0.95, 1.04]	-0.01	[-0.08, 0.07]	0.04	[-0.11, 0.19]	0.00	[-0.12, 0.11]
is_B:pushing_self_cw	0.92	[ 0.79, 1.09]	1.02	[ 0.96, 1.08]	-0.06	[-0.14, 0.03]	-0.07	[-0.25, 0.10]	-0.05	[-0.18, 0.08]
is_A:pushing_partner_cw	1.03	[ 0.89, 1.20]	1.03	[ 0.97, 1.09]	0.00	[-0.08, 0.09]	0.00	[-0.18, 0.17]	0.00	[-0.13, 0.13]
is_B:pushing_partner_cw	1.00	[ 0.87, 1.15]	0.98	[ 0.93, 1.02]	0.05	[-0.02, 0.13]	-0.09	[-0.25, 0.06]	0.05	[-0.09, 0.18]
is_A:day	0.73	[ 0.48, 1.12]	1.08	[ 0.91, 1.28]	-0.06	[-0.33, 0.21]	0.06	[-0.51, 0.64]	0.30	[-0.24, 0.85]
is_B:day	1.07	[ 0.69, 1.65]	0.94	[ 0.78, 1.12]	-0.01	[-0.26, 0.25]	-0.55	[-1.13, 0.05]	0.03	[-0.58, 0.62]
is_A:weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cw	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cw	1.03	[ 0.93, 1.14]	1.01	[ 0.98, 1.04]	-0.02	[-0.07, 0.04]	0.01	[-0.10, 0.12]	0.08	[-0.03, 0.19]
is_B:barriers_self_cw	0.99	[ 0.91, 1.08]	0.98	[ 0.95, 1.00]	0.02	[-0.03, 0.06]	-0.03	[-0.13, 0.06]	0.01	[-0.10, 0.11]
Between-Person Effects
is_A:pushing_self_cb	0.47	[ 0.19, 1.15]	0.81	[ 0.56, 1.19]	0.42	[-0.22, 1.08]	1.01	[-0.47, 2.54]	-0.40	[-1.15, 0.37]
is_B:pushing_self_cb	0.73	[ 0.35, 1.52]	0.77	[ 0.52, 1.13]	-0.28	[-0.84, 0.26]	-0.80	[-1.96, 0.41]	-0.76	[-1.52, 0.04]
is_A:pushing_partner_cb	1.31	[ 0.68, 2.60]	1.12	[ 0.77, 1.64]	-0.07	[-0.61, 0.49]	-1.21	[-2.51, 0.05]	0.10	[-0.63, 0.86]
is_B:pushing_partner_cb	2.03	[ 0.84, 4.97]	0.74	[ 0.48, 1.16]	0.41	[-0.23, 1.06]	1.16	[-0.36, 2.58]	0.51	[-0.54, 1.59]
is_A:weartime_self_cb	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_B:weartime_self_cb	NA	NA	1.00	[ 1.00, 1.00]	NA	NA	NA	NA	NA	NA
is_A:barriers_self_cb	1.08	[ 0.89, 1.31]	1.07	[ 0.99, 1.17]	0.01	[-0.14, 0.16]	0.21	[-0.11, 0.54]	0.10	[-0.11, 0.32]
is_B:barriers_self_cb	1.04	[ 0.82, 1.30]	0.98	[ 0.89, 1.09]	-0.07	[-0.25, 0.11]	0.04	[-0.35, 0.42]	-0.16	[-0.46, 0.14]
Random Effects
sd(is_A)	0.63	[ 0.39, 0.97]	0.34	[0.24, 0.47]	0.79	[0.61, 1.02]	1.43	[1.08, 1.88]	0.16	[0.01, 0.42]
sd(is_B)	0.47	[ 0.29, 0.68]	0.34	[0.23, 0.48]	0.68	[0.51, 0.89]	1.52	[1.15, 1.98]	0.26	[0.01, 0.64]
sd(is_A:pushing_self_cw)	0.12	[ 0.01, 0.33]	0.09	[0.01, 0.18]	0.07	[0.00, 0.18]	0.15	[0.01, 0.37]	0.29	[0.11, 0.48]
sd(is_B:pushing_self_cw)	0.46	[ 0.14, 0.85]	0.11	[0.01, 0.20]	0.12	[0.01, 0.27]	0.19	[0.01, 0.49]	0.23	[0.04, 0.43]
sd(is_A:pushing_partner_cw)	0.43	[ 0.11, 0.84]	0.11	[0.02, 0.20]	0.17	[0.04, 0.32]	0.28	[0.02, 0.60]	0.13	[0.00, 0.39]
sd(is_B:pushing_partner_cw)	0.11	[ 0.01, 0.30]	0.03	[0.00, 0.09]	0.10	[0.01, 0.24]	0.11	[0.00, 0.31]	0.12	[0.01, 0.35]
Additional Parameters
ar[1]	0.04	[-0.92, 0.94]	0.21	[0.14, 0.28]	0.18	[0.12, 0.24]	0.24	[0.18, 0.30]	0.13	[0.00, 0.27]
nu	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
shape	0.43	[ 0.40, 0.47]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sderr	0.05	[ 0.00, 0.13]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]	NA	[ NA, NA]
sigma	NA	[ NA, NA]	0.52	[0.50, 0.54]	0.85	[0.82, 0.89]	1.80	[1.73, 1.88]	1.00	[0.92, 1.08]

# prepare export to excel
all_models$Pushing_apim <- prepare_for_export(pushing_all_apim_models)

openxlsx::write.xlsx(all_models, file = "Output/AsPreregistred/AllModels.xlsx")